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ABSTRACT 

The  free  boundary  problem  for  a fully  penetrating  well  of  radius  r 

and  filled  with  water  to  a depth  h,  in  a layer  of  soil  of  depth  H,  radius 

R and  permeability  k(x,y)  can  be  formulated  as  follows:  Find  < pe  C^[r,R]  and 
2 — 

u c C (ft)  n C (ft)  such  that  (xku  ) + (xku  ) = 0 in  ft,  u(R,y)  = H for 

x x y y 

0<y<H,  u (x,0)  = 0 for  r<x<R,  u(r,y)  = h for  0<y<h,  u(r,y)  = y 
n 

for  h<y<  ip  ( r)  , u (x,^  (x) ) = 0 and  u(x,  v>  (x) ) = 'P  (x)  for  r < x < R,  where 

n 

ft  = { (x,y)  : r<x<R,  0<y<  >Mx) } . The  res\ilts  of  Benci  [Annali  di  Mat. 

100  (1974) , 191-209]  are  used  to  derive  a variational  inequality  and  to  prove 
existence  and  uniqueness.  The  problem  is  approximated  using  piecewise  linear 
finite  elements  and  0(h)  convergence  of  the  approximate  solutions  is  proved 
using  recent'  results  due  to  Brezzi,  Hager,  and  Raviart. 
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EXPLANATION 

When  an  axisymmetric  well  is  sunk  in  a water  aquifer  the  water  flows 
through  the  ground  towards  the  well  where  it  can  be  removed  by  pumping. 

After  a sufficiently  long  time  the  flow  becomes  steady.  We  analyse  this 
problem  using  the  assumptions  of  saturated  - unsaturated  porous  flow: 
the  upper  part  of  the  aquifer  is  denuded  of  water  and  is  dry;  the  lower  part 
of  the  aquifer  is  saturated  (wet)  and  in  it  the  flow  is  governed  by  a 
linear  second  order  elliptic  differential  equation  which  is  derived  from  the 
equation  of  continuity  and  Darcy’s  law  (a  linear  relationship  between  the 
rate  of  flow  and  the  water  pressure  gradient) . The  major  mathematical 
difficulty  is  that  the  interface  between  the  dry  and  wet  regions  is 
unknown  - it  is  called  a free  boundary.  We  reformulate  the  problem  as  a 
variational  inequality,  that  is,  a variational  problem  in  which  the  solution 
is  subject  to  inequality  constraints  (in  our  case,  the  solution  must  be 
non-negative) . We  prove  existence  and  uniqueness  (which  were  previously 
unknown) , introduce  a numerical  scheme  based  on  finite  elements,  and  prove 
convergence  of  the  numerical  approximation  to  the  exact  solution.  Numerical 
results  are  given  for  the  case  of  constant  permeability  (homogeneous  isotropic 
soil)  and  the  computer  program  is  listed. 
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THE  NUMERICAL  SOLUTION  OF  AXISYMMETRIC  FREE 
BOUNDARY  POROUS  FLOW  WELL  PROBLEMS  USING 
VARIATIONAL  INEQUALITIES 

Colin  W.  Cryor  and  Hans  Fetter 

1.  Introduction 

The  steady  state  problem  to  be  considered  is  shown  in  Figure  1.1.  An  axisymmetric 
well  of  radius  r is  sunk  into  a layer  of  soil  of  depth  H and  radius  R . The  bottom 
of  the  soil  layer  is  impervious.  The  outer  boundary  of  the  soil  adjoins  a catchment  area 
and  the  hydraulic  head  u is  equal  to  the  constant  H along  this  boundary.  The  water 
seeps  towards  the  well  and  a pump  (not  shown)  maintains  the  water  level  in  the  well  at 
a constant  height  h^  . The  water-air  interface  is  a free  boundary  which  intersects  the 
well  wall  at  a height  hs  . 

The  mathematical  problem  can  now  be  formulated  as  follows  (see  Hantush  [1964]  , 

Bear  [1972],  and  Cryer  [1976,  p.  86]): 

Problem  A (Classical) 

Find  functions  V- (x)  (the  height  of  the  free  boundary)  and  u(x,y)  (the  hydraulic 
head)  such  that  (from  the  equation  of  continuity  and  Darcy's  law): 

div  (k  grad  u)  = > + > = °>  in  , (1.1) 

together  with  the  boundary  conditions, 


u = H, 

on 

ab (r^> , 

(constant  hydraulic  head)  , 

(1.2) 

3u 

3n  “ °' 

on 

bc(t.)  , 

4 

(impervious  boundary)  , 

(1.3) 

u = h , 

w 

on 

cD(r2>, 

(interface  with  water  at  rest) , 

(1.4) 

u = y. 

on 

DE(I^), 

(interface  with  air)  , 

(1.5) 

u = y. 

on 

EA(r0), 

(interface  with  air)  , 

(1.6) 

£ 

, on  EA(rQ), 

(streamline) . 

(1.7) 
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Here,  0 is  the  (unknown)  domain, 


t 


fl  " { (x,y)  : 0 < y < ip  (x)  , r<x<R  } , 

£ 

and  denotes  the  outward  normal  derivative.  Finally,  k = k(x,y)  = xx  where  the 

permeability  of  the  soil  is  denoted  by  K = K(x,y)  . It  is  assumed,  that  k is  of 
the  form 


k(x,y)  = exp[f(x)  + g(y)] 

where  f Cx)  and  g(y)  are  continuously  differentiable  and 

g'  (y)  >_0  . D 

In  particular  if  the  permeability  K is  constant,  K = 1 say,  then 

k(x,y)  = x = exp[£,n  x] 

so  that 

f(x)  = Inx;  g (y)  =0  . 

We  will  later  use  the  fact  that 

u = y + p/pg  , 


(1.8) 


(1.9) 


(1.10) 


(1.11) 


where  g is  the  acceleration  due  to  gravity,  p is  the  fluid  pressure,  and  p is  the 
fluid  density. 

Ever  since  C.  Baiocchi  [1971]  introduced  a mathematically  rigorous  and  from 
the  numerical  point  of  view  efficient  approach  to  the  solution  of  various  free  boundary 

problems  related  to  fluid  flow  through  porous  media,  numerous  studies  have  appeared 
in  the  literature  which  extend  his  results  in  many  directions;  Cryer  [1 976]  and  Baiocchi, 
Brezzi,  and  Comincioli  [1976]  give  bibliographies. 

The  basic  idea  introduced  by  Baiocchi,  can  be  summarized  as  follows:  Through  a 
suitable  change  of  the  unknown  variable,  the  free  boundary  problem  is  reduced  to  that  of 
minimizing  a quadratic  functional  on  a closed  convex  set.  This  reformulation  of  the  problem 
not  only  enables  one  to  determine  various  properties  of  the  solution,  but  it  also  offers 
the  advantage  that  the  new  problem  can  readily  be  solved  numerically  by  several  methods 
including  finite  differences  and  finite  elements. 
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The  first  problem  considered  by  Baiocchi  was  the  "model  problem"  of  porous  flow 


between  two  water  reservoirs  of  different  levels  separated  by  a homogeneous  rectangular 
dam  with  horizontal  impervious  base.  Subsequently,  Baiocchi,  Comincioli,  Guerri , and  Volpi 
(1973)  and  Baiocchi,  Comincioli,  Magenes,  and  Pozzi  [1973]  considered  the  case  of  a rec- 
tangular dam  consisting  of  two  homogeneous  layers,  either  horizontal  or  vertical.  A 
further  development  is  due  to  Benci  [1973,  1974]  who  considered  a rectangular  dam  with 
a permeability  coefficient  of  the  form 

tc(x , y ) = exp(f(x)  + g (y ) ) . (1.12) 

The  problem  of  a fully  penetrating  axially  symmetric  well  (Figure  1)  s the  axially 
symmetric  equivalent  of  the  problem  of  flow  through  a rectangular  porous  dam.  The  well 
problem  is  of  considerable  technological  importance.  Hantush  [1964]  gives  a lengthy 
survey,  and  Cryer  [1976,  p.  86]  gives  further  references. 

Polubarinova-Kochina  [1962,  p.  283],  Hantush  [1964,  p.  362],  and  Bear  [1972, 
p.  368]  show  that  for  constant  permeability  the  rate  of  flow  is  given  by 

Q w = HK:(H2-h2)/£n(R/r)  . (1.13) 

Mauersberger  [1963]  and  Youngs  11971]  obtain  exact  expressions  in  closed  form  for  the 
rate  of  flow  with  variable  permeability  of  the  form  (1.12).  Mauersberger  [1965]  has 
derived  a variational  principle.  However,  apart  from  the  not  entirely  rigorous  results 
of  Mauersberger  [1965a] , there  appear  to  be  no  results  on  existence  and  uniqueness.  In 
this  paper  the  results  of  Benci  [1973,  1974]  are  applied  to  the  well  problem  to  obtain 
existence  and  uniqueness.  A finite  element  method  is  then  used  to  obtain  numerical  ap- 
proximations, and  error  estimates  are  obtained. 

After  completing  this  report  we  learned  that  the  case  of  constant  K had  been 
formulated  as  a variational  inequality  and  solved  numerically  by  Elliott  [1976,  p.  62]  . 


2.  Notation  and  preliminaries 

We  state  here  certain  basic  definitions  and  results  about  Sobolev  spaces.  This 

material  can  be  found  in  Adams  11975], 

2 

If  ft  is  an  open  set  in  R the  following  spaces  may  be  constructed: 

C(ft)  = C°(ft)  is  the  Banach  space  of  functions  which  are  bounded  and  uniformly  continuous 
on  ft  with  norm 

||  u;  C (ft)  ||  = max  | u | . 

ft 

Cm(ft)  is  the  vector  space  of  functions  u which,  together  with  all  their  partial  derivatives 
D°u  of  order  |a|  < ra  are  continuous  on  ft  . Here  a = (a, , . . . ,a  ) with  a.  a 
non-negative  integer, 

M - l Kl  . 

j=i  3 

and 


a 

D u 


3'a'u 


a a 

8x  . • . 3x  n 
1 n 


cm(ft) 


anach  space  of  those  functions  u CC  (ft)  for  which  Du  is  bounded 
and  uniformly  continuous  for  0 <Ja|  m . If  u f C (ft)  then  D u can  be  ex- 

tended continuously  to  ft  and  we  may  regard  Du  as  being  defined  on  ft  , and 
it  is  readily  shown  that  Dau  is  bounded  and  uniformly  continuous  on  ft  . 

Cm(ft)  is  equipped  with  the  norm 


|u;Cm(ft) 


max  sup  |Dau(x)|  , 
0<  |a|£m  xeft 


max  max  |Dau(x)| 
0<Ja|<^  m x c ft 


C (ft)  is  the  vector  space  consisting  of  functions  u such  that  u*  C (ft)  for  all  m , 
0 < m < “ . 
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c™(ft)  is  the  subspace  of.  cm (H)  consisting  of  functions  u which  "vanish  near  8ft"  , 
that  is,  the  functions  ue  cm(ft)  with  compact  support  in  ft  . C™(ft)  is  of  course 
also  a subspace  of  , 

is  the  subspace  of  C°° (ft)  consisting  of  functions  ueC  (ft)  with  compact  support 
in  ft  . 

LP(ft)  is  the  Banach  space  of  real  measurable  functions  u defined  on  ft  for  which 


llul^  H || u;  (ft)  ||  = [fn  |u(x)  |Pdx]1/p  < “ • 

Here,  l<p<  Two  elements  of  LP(ft)  are  identified  if  they  are  equal  a.e. 

(almost  everywhere) . 

Hm,P(ft)  is  the  completion  of  {uf  cm(fi):  ||u  ||  ^ < °°)  with  respect  to  the  norm 


( l (||Dau||  )P)VP  . 

0<  |a|_<  m 1 


Hm,P  ft  is  the  closure  of  C (ft)  , in  II  (ft)  . 

0 0 

Hm(ft)  = Hm'2(ft)  . The  norm  in  H™  is  denoted  by  ||  • II  m • 

H™(ft)  = Il™'2(ft)  . By  the  generalization  of  Poincare's  Lemma,  if  ft  is  bounder  then  the 


norm 


m,p 


H™,P(ft) 


is  equivalent  to  the  norm 


defined  bv 

m,p 


■m,p 


= I1 

a =m 


|Dau 


)P]Vp 


The  spaces  Hm,P(ft)  and  H™,P(ft)  are  called  Sobolev  spaces  or  Beppo-Levi  spaces. 
There  is  a useful  alternative  definition  of  the  Sobolev  spaces.  If  uf  L (ft)  and 
if  for  some  a there  exists  v^  e L^(ft)  such  that 

/ u(x)Da<Mx)dx  = (-1)^  / v (x)  <|)  (x)  dx  , 

ft  « 

for  all  $ c c“(ft)  , 
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Weak  formulation  of  the  problem 


I 


In  the  (classical)  Problem  A it  is  implicitly  assumed  that  the  free  boundary  fo 
and  solution  u are  "sufficiently  smooth".  Here  we  reformulate  the  problem  as  a second 
problem,  Problem  B,  in  which  these  smoothness  assumptions  are  relaxed. 

2 — 

Let  ta.u}  be  a "smooth"  solution  of  Problem  A.  We  cannot  assume  that  ut  C (ft) 

because  -r—  is  discontinuous  at  the  corner  D (see  Figure  1.1).  It  might  also  appear 
3y 

that  u could  be  discontinuous  at  the  corners  B and  C;  this  is  not  the  case  since  by 
reflecting  both  u and  ft  in  the  x-axis  we  obtain  a solution,  u*  say,  in  the  reflected 
domain,  ft'  say.  Let 

ft"  = ft  uH'tiT,  , 

4 


and  let 


u" 


u,  on  ft  u 

u 1 , on  ft 1 


Then  3(2"  is  smooth  at  B and  C so  that  u"  is  well-behaved  at  B and  C . We  thus 

assume  that  "sufficiently  smooth"  implies  the  following:  the  free  boundary  has  a 

2 — 

continuous  outward  normal  n ; u e C (ft)  n C(ft);  and  u is  continuously  differentiable 
in  a neighborhood  of  u T . 

To  obtain  the  weak  formulation  of  the  problem  we  proceed  as  follows.  For  any 
ip  e C (ft)  which  vanishes  in  a neighborhood  of  1’  U u it  follows  from  Green  s 

theorem  (see  Appendix  C)  and  (1.1)  that 

/ k grad  u grad  tj)  dxdy 
ft 

= f \)j  k |—  ds  - / i|i  div(k  grad  u)dxdy  , (3.1) 

3ft  ft 

= 0 . 
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Let 


V = (vf  H1  (ft) : v = 0 on  Tj  u T2  u r j}  ; (3.2) 

that  is,  V is  the  closure  in  (SI)  of  the  functions  i e C (ft)  which  vanish  in  a 

neighborhood  of  u I"2  u rj  . Since  every  v e V is  the  limit  in  (fl)  of 

functions  t|)  for  which  (3.1)  holds,  we  have  that 

/ ^ k grad  u grad  v dxdy  = 0,  for  all  veV  . (3.3) 

Thus,  { V>,  u)  is  a solution  of 


Problem  B (Weak) 

Find  <f  e C[r,Rj  and 
boundary  conditions 


u e H*(ft)  n C(ft)  such  that 

u = H,  on  , 

u = h , on  r„  , 
w 2 

u = y , on  , 

u = y,  on  rQ  . □ 


u satisfies  (3.3)  and  the 

(3.4) 

(3.5) 

(3.6) 

(3.7) 


Formulation  as  a variational  inooualit' 


In  this  section  we  follow  Denci  11973,  1974]  and  reformulate  the  weak  problem, 


Problem  B,  as  a variational  inequality,  Problem  C.  The  derivation  given  here  is  in- 


tended to  complement  that  of  Benci : it  is  hoped  that  the  derivation  given  here  brings 


out  more  clearly  the  essential  steps  in  the  derivation. 


The  first  step  is  to  introduce  a "Baiocchi  function"  w(x,y)  defined  on  the 


rectangle 


D = { (x,y)  : r£x£R,  0 < y < h)  , 


as  follows: 


w(x,y)  = 


/ exp(g(t))  fu(x,t)  - t]dt,  for  (x,y)  e fi  , 


0,  for  (x,y)  c D - ft  . 


We  show  that  w satisfies  a differential  equation  at  points  (x,y)  e f2  . Dif- 


ferentiating (4.2)  with  respect  to  x we  find  that 


w (x,y)  - / exp[g(t)]  u (x,t)dt  + 


+ ^’(x)  exp[g  (*?  (x) ) ] [u(x,v>(x))  — (x)  J , 


<t  (x) 

/ exp(g(t)]  u (x,t)dt.  , 

y 


since,  by  (1.6),  u(x,  (x)  ) - v>  (x) 


Multiplying  (4.3)  by  exp[f(x)]  and  then  differentiating  with  respect  to  x we  obtain 


(exp(f  (x)  ] wx<x,y)) 


/ [exp [ f (x)  + g (t)  ] u (x,  t)  ] dt  + 

y xx 


+ 'fi ' (x)  exp[f(x)  +g  («p  (x) ) J u (x,  ?(x)) 


We  now  carry  out  similar  computations  for  the  derivatives  of  w with  respect 


to  y . 


I 


[ 

[ 

i 


Wy(x,y)  = -explg(y) J |u(x,y)  - y]  . (4.5) 

Multiplying  by  expl-g(y)]  and  differentiating  with  respect  to  y we  obtain, 

(expl-g(y)]  w^(x,y))^ 

= 1 - uy(x,y)  , 

= 1 - exp  (-g  (y)  ) (exp[g(y))  uy(x,y)]  , 

V>  (x) 

= 1 + cxp[-g  (y)  ] {/  (exp[g(t)  ] ut(x,t))t  dt  - 

y 

- explgte (x) ) ] Uy(x,  <p  (x)  ) } . (4.6) 

Multiplying  (4.4)  by  exp[-g(y))  and  (4.6)  by  exp[f(x))  and  adding,  we  obtain  that 

(exp[f  (x)  - g(y))  w (x,y))  + (exp[f(x)  - g(y))  w (x,y)) 

x x y y 

= exp[f (x) ] + 

V>  (x) 

+ exp[-g  (y)  ] / [(k(x,t)  u^txjt))^  + (k(x,t)  ut<x,t))t)dt  + 

y 

+ exp[-g(y)  + f (x)  + g(tp(x))]  lv>'(x)  u (x,v’(x))-  u (x,^(x)))  . (4.7) 

x y 

Since  u satisfies  div(k  grad  u)  = 0,  and  since  on  the  free  boundary 

2 1/2 

V>'  u - u = -[1  + (*')]'  u = 0 , 

x y n 

the  last  two  terms  on  the  right  hand  side  of  (4.7)  vanish.  Thus,  if  the  elliptic  operator 
L is  defined  by 

(Lw)  (x,y)  = (exp[f (x)  - g(y)]w  ) + (exp[f(x)  - g(y)]w  ) , (4.8) 

xx  y y 

we  have  from  (4.7)  that 

Lw  = exp[f(x)],  in  U . (4.9) 

Since,  from  (4.2),  w = 0 in  D - U,  we  also  have  that 

Lw  = 0,  in  D - fl  . (4.10) 

At  this  point  it  is  appropriate  to  observe  that,  from  (4.2),  (4.3),  and  (4.5)  , 
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(4.11) 


w -w  ~ w = 0,  on  r„  . 
x y n 0 

We  denote  by  $ the  boundary  values  of  w . The  function  $ can  be  determined  as 
follows.  From  (1.2)  and  (4.2)  it  follows  immediately  that 

<Mx,H)  = 0,  on  I’g  . (4.12) 

H 

4>(R,y)  = / exp[g(t)]  (H-tldt,  on  r . (4.13) 

y 

From  (1.4)  and  (1.5)  we  see  that 

h 

w 

<Mr,y)  = / exp[g(t)  ] [h  -t]dt,  on  T.  , (4.14) 

v/  z 

y 

4>(r,y)  = 0,  on  F^  . (4.15) 

To  determine  $ on  we  must  proceed  more  indirectly.  From  (4.4)  we  see  that 

Iexp[f(x)]  w (x , 0) ) 
xx 

<e  (x) 

= / Iexp[f(x)  + g(L)]  u (x,t)]  dt  + 

0 xx 

+ ' (x)  cxp[f(x)  + g (sP  (x) ) 1 u^fx.ipfx))  . 

Replacing  the  integrand  using  the  equation  div  k grad  u = 0 and  then  integrating,  wc 
obtain 

[exp[f (x)  WX(X<°)HX 
<P  (x) 

= -/  [exptf(x)  + g (t) ) u (x,t)J  dt  + 

0 t 

+ V>'(x)  exp[f(x)  + g (v»  (x) ) } u^fx.vMx))  , 

= exp(f (x)  + g(0 ))  u^(x,0)  + 

+ exp[f  (x)  + g(v>(x))]  N>*  (x)  u (x,i p (x) ) - u (x,tf(x))]  , 

x y 

= 0 , 

since,  by  (1.3)  and  (1.7),  u^  = 0 on  F^  and  u^  = 0 on  F^  . 
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We  thus  conclude  that  on  f4  $ satisfies  the  linear  second  order  differential  equation 


(explf  (x)  ]♦  ) = 0 > 


(4.16) 


with  general  solution 


4>(x,0)  = A + B / exp|-f(t)]dt 


(4.17) 


Since  the  value  of  $(x,0)  is  given  at  x = r by  (4.14)  and  at  x = R by  (4.13)  , the 

values  of  A and  B can  be  determined.  We  obtain: 


4>(x,0)  = 4><r,0)  + (i>(R,0)  - i>(r,0)]  F(x)/F(R)f  on  r , 


(4.18) 


Fix)  = / exp  I - f ( t.) ) dt 


(4.19) 


From  (1.11), 


u = y + p/pg 


(4.20) 


where  g is  the  acceleration  due  to  gravity,  p is  the  fluid  pressure,  and  p is  the  fluid 
density.  The  fluid  pressure  p is  non-negative  for  physical  reasons  so  that  u-y2_0 
and  hence  w>0  . Unfortunately,  this  cannot  be  proved  mathematically  without  some 
additional  assumptions.  One  approach  is  as  follows.  Consider  the  function  v = y - u . 


div(k  grad  v)  = div  k grad  y - div  k grad  u 


(k)  , 


= explf  (x)  + g(y)l  g' (y) 


div  k grad  v > 0 in  fJ 


provided  that  g'  (y)  > 0 . From  the  maximum  principle  for  elliptic  equations  (Courant 
and  Hilbert  (1962,  p.  326))  we  conclude  that  v attains  its  maximum  in  fl  on  30  • 
However,  from  the  boundary  conditions,  v £ 0 on  Tp  U ^ u u ^ . Also,  on 
we  have  that  -13- 


V = 1 - u =1 

y y 

so  that  v cannot  attain  its  maximum  on  . Consequently,  we  see  that  v 0 in 


ft  . That  is,  wc  have  shown  that 


provided  that 


> 0 in  ft 


g*  (y)  > 0 . 


Let  a be  the  bilinear  operator  defined  on  H3  (D)  X H3(D)  by, 


i(u,v)  = / exp[f(x)  - g(y)]grad  u grad  v dxdy  , 
D 


5 / exp[f(x)  - g (■/))  (u  v + u v Jdxdy 


xx  y y 


Let  j be  the  linear  functional  defined  on  H (D)  by, 

j (v)  = / exp[f(x)]v  dx  dy 
D 


(4.21) 


(4.22) 


(4.23) 


(4.24) 


Let  K be  the  closed  convex  set 


K = {v  e II3  (D)  : v - w c H3  (D)  and  v 0 a.e.  in  D) 


(4.25) 


Now  let  v be  any  element  in  K . Then 


a(w,v-w)  + j(v-w) 


f exp[f(x)  - g(y)]  grad  w grad  (v-w)dxdy  + 

D 


+ J exp|f(x))  (v-w)dxdy  , 

D 

/ exp[f(x)  - g(y)l  grad  w grad  (v-w)dxdy  + 
ft 

+ j exp[f(x)]  (v-vj)dxdy  , 

ft 

since  w = 0 in  D - ft  . Integrating  by  parts,  and  noting  that  either  v-w  =0  or 

* 0 on  3ft  , we  obtain  that 
3n 


I 


a(w,v-w)  + j (v-w) 


/ I-Lw  + explf(x)]]  (v-w)dxdy, 

n 


from  (4.9)  . 


We  have  thus  shown  that  if  g' (y)  > 0 then  w satisfies  the  variational  in- 


equality 


a(w,v-w)  + j (v-w)  > 0 


(4.26) 


More  precisely,  Benci  973,  1974)  proves 


Theorem  4 . 1 


If  u is  a solution  of  the  weak  problem  and  g' (y)  ^ 0 then  w c H (D)  n C(D) 


and  w satisfies  the  variational  inequality:  Find  w e K such  that 

a(w,v-w)  + j (v-w)  > 0 , 


(4.27) 


for  all  v r K . 


Remark  1 


Benci  11974,  p.  194]  and,  in  a special  case,  Baiocchi,  Comincioli,  Magenes,  and 
Pozzi  [1973,  p.  19)  assert  that  if  u is  a solution  of  the  weak  problem  and 


u-y,  in  fi  , 


0 , in  D-0  , 


then  IT  e H (D)  . We  have  been  unable  to  "follow  their  arguments,  and  it  seems  to  us 
that  their  arguments  require  that  0 be  such  that  trace  theorems  can  be  applied.  This 
is  not  a serious  difficulty,  since  if  u is  "sufficiently  smooth"  then  we  certainly 
have  that  n t (D)  . 

Remark  2 

Both  Benci  (1974)  and  Baiocchi,  Comincioli,  Magenes  and  Pozzi  [1973)  assume  that 
t is  a monotonically  decreasing  function.  While  undoubtedly  true,  this  assump- 
tion does  not  seem  to  be  necessary  for  the  derivation  of  the  variational  inequality. 
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Remark  3 


We  have  assumed  that 
tion  can  be  relaxed  somewhat: 


f and  g are  continuously 
Benci  [1974]  only  requires 


f c H1'2^  [r.Rj  . 


g e H1,2+U  [0,H]  , 


differentiable. 

that 


This  assump- 


for  some  p > 0 . 

The  assumption  that  f be  reasonably  smooth  is  not  unduly  restrictive.  However, 
the  soil  around  a well  often  consists  of  N horizontal  layers  of  different  constant 
permeabilities,  and  g is  then  not  in  [0,H]  . 

Rema  rk  4 

While  condition  (4.22)  is  not  necessary,  it  appears  that  some  condition  must  be 
imposed  upon  g (y ) . Baiocchi  and  Friedman  [to  appear]  refer  to  numerical  solutions  of 
the  variational  inequality  (4.27)  by  Comincioli  showing  that  some  choices  of  g ap- 
parently lead  to  negative  water  pressures  p,  and  this  has  recently  been  proved  by 
Friedman  [1977J . 


Remark  5 

For  results  when  k is  not  of  the  form 


see  Baiocchi  [1976]. 


k = exp[f  (x)  + g(y)  ] 


Remark  6 

Although  it  is  not  apparent  from  the  above  derivation  of  the  variational  inequality 
(4.27),  the  fact  that  the  boundary  values  4>  of  the  Baiocchi  function  w can  be  ex- 
plicitly computed  is  related  to  the  fact  that,  as  shown  by  Mauersberger  [1969]  and 
Youngs  [1971] , the  flow  rate  Qw  can  bo  explicitly  computed. 

If  a problem  is  such  that  cannot  be  determined  explicitly,  then  an  additional 

complication  is  introduced,  namely  that  the  unknown  Q must  also  be  found;  see 
Baiocchi,  Comincioli,  Magenes,  and  Pozzi  [1973,  p.  46]. 
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5.  Properties  of  the  solution 

By  assumption  f (x)  and  g(y)  are  bounded.  Thus,  the  bilinear  functional  a and 
the  linear  functional  j are  continuous  on  H* (D) : 


a(v,v)  < ct2<||  v||  l 2)‘ 


h(v>|  < e2||  v||  i 2 


(5.1) 

(5.2) 


for  all  v e H1(D).  where  c*2  and  $2  are  constants. 

Remembering  that  the  norms  | • |m  p and  ||  • ||  ^ are  equivalent  on  H™'P(D)  (see 
section  2)  we  can  conclude  that  a is  coercive  on  h*(d): 


a(v,v)  > ( ||  v ||i2)  , 


(5.3) 


for  all  v £ H*(D),  where  is  a strictly  positive  constant. 

It  follows  from  the  basic  theory  of  variational  inequalities  (Stampacchia  [1964])  that 
Theorem  5. 1 

There  exists  a unique  solution  w £ H*(D)  of  the  variational  inequality  formulation 
(4.27)  of  the  axisymmetric  well  problem.  D 

Although  Theorem  5.1  answers  the  most  basic  questions,  namely  regarding  existence 
and  uniqueness,  there  remain  a number  of  interesting  questions  which  we  now  mention: 

1.  How  smooth  is  w ? 

If  w is  to  be  a solution  of  the  classical  problem,  Problem  A,  then  w cannot 
just  be  in  H (r»  . Moreover,  as  will  be  seen  in  section  6,  the  smoothness  of  w plays 
an  important  role  in  the  error  analysis  of  approximation  methods. 

It  follows  from  the  results  of  Benci  [1974,  p.  200)  that 


w c H2,p(d)  , 


(5.4) 


for  any  p satisfying  1 < p < <®,  since  D has  the  cone  property  it  is  a consequence 
of  the  Sobolev  embedding  theorem  (Adams  [1975,  p.  97))  that 

w £ C1 (D)  . (5.5) 
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r 


p 


* 


It  would  be  of  considerable  interest  for  numerical  applications  if  it  could  be 

shown  that  w was  even  smoother.  In  particular,  if  w c II^2  c,p  for  any  c > 0 then 
the  quadratic  approximations  due  to  Brezzi,  Hager,  and  Raviart  [to  appear]  could  be  ap- 
plied. At  the  time  of  writing  we  do  not  know  whether  this  is  true.  Of  course,  we  can- 
not expect  w to  be  very  smooth  because  w has  a discontinuity  across  the  free  boundary. 
References  on  the  regularity  of  solutions  of  variational  inequalities  include:  Brezis 
[1971];  Brezis  and  Kinderlehrer  [1973/74];  Brezis  and  Stampacchia  [1968];  Lewy  and 
Stampacchia  [1969]. 

2.  What  are  the  properties  of  w and  0 ? 

Many  interesting  questions  suggest  themselves  concerning  w and  [2.  In  particular 
the  following  properties  are  known: 

(a)  w <0  and  w < 0 in  D . (Benci  [1974,  p.  200  and  p.  202]) 

x — y — 

(b)  is  continuous  and  strictly  decreasing.  (Benci  [1974,  p.  207]  and  Baiocchi 
and  Friedman  [to  appear]  .) 

(c)  For  other  results  for  the  case  k = 1 see  Caffarelli.  [to  appear]  and  Jensen 
[to  appear] . 
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6.  Numerical  approximation 


It  has  been  shown  in  the  previous  sections  that  the  Baiocchi  function  w satisfies 
the  variational  inequality  (4.27):  Find  w t K such  that  for  all  v c K , 

a(w,v-w)  + j (v-w)  0 . (6.1) 

Since  a is  symmetric,  that  is 

a(v,w)  = a(w,v),  for  all  v.w  e V , 

there  is  a connection  between  the  variational  inequality  (6.1)  and  the  unilateral  mini- 
mization problem 


Min  J(v)  , 

V£K  (6.2) 

J(v)  = a (v,  v)  + 2j  (v)  . 

This  connection  is  given  by  the  following  theorem  which  is  well-known  but  which  we  prove 
for  the  convenience  of  the  reader. 

Theorem  6.1 

Let  a(v,w)  be  a symmetric  bilinear  form  satisfying  a(v,v)  > 0 for  all  v e V . 
Then  w is  a solution  of  the  variational  inequality  (6.1)  iff  w is  a solution  of  the 
unilateral  minimization  problem  (6.2). 

Proof.  Suppose  that  u is  a solution  of  (6.1).  Then,  for  any  v e K , 

J (v)  = a (v,v)  + 2j  (v)  , 

= a(u+(v-u),  u + (v-u))  + 2 j (u  + (v-u) ) , 

= J(u)  + 2[a(u,v-u)  + j(v-u)]  + a (v-u,  v-u)  , 

_>  J(u) 

so  that  u is  a solution  of  (6.2). 

Now  suppose  that  u is  a solution  of  (6.2).  Then,  since  K is  convex,  for  any 
v e K we  have 

u + t (v-u)  e K,  for  0 t 1 . 

Thus 
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G (t)  = J(u+t(v-u))  - J(u,u) 


2t[a(u,v-u)  + j(v-u)]  + t a(v-u,v-u) 


>;  0,  for  0 < t < 1 . 

It  follows  that 

G' (0)  = 2[a(u,v-u)  + j(v-u)]  >_  0 ; 
so  that  u is  a solution  of  (6.1).  D 

We  note  that  Theorem  6.1  does  not  assert  that  either  the  variational  inequality  or 
the  unilateral  minimization  problem  has  a solution.  In  the  present  case  the  bilinear 
form  is  coercive  (see  (6.3))  and  wo  know  from  Theorem  5.1  that  a solution  exists  and  is 
unique.  (Theorem  6.1  with  the  added  assumption  of  coercivity  is  given  by  Lions  [1971, 
p.  9]). 

We  approximate  w by  choosing  a finite-dimensional  approximation  and  solving 

the  finite-dimensional  problem:  Find  w^  c 

J(wh>  = Min  J (v ^)  . (6.3) 


The  convex  set  is  constructed  as  follows.  The  domain  D is  triangulated  as 

shown  in  Figure  6.1. 


Figure  6.1:  The  triangulation  of  D 
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The  subdivisions  are  not  necessarily  uniform,  but  it  is  assumed  that  there  is  a constant 
6 > 0 such  that 

(maximum  interval  lengtH  < h < 3 (minimum  interval  length)  < (6.4) 

p - — 

where  h is  a measure  of  the  fineness  of  the  subdivision.  The  set  of  interior  gridpoints 

will  be  denoted  by  D and  the  set  of  boundary  gridpoints  will  be  denoted  by  3d  . 

h h 

We  denote  by  V.  the  space  of  piecewise  linear  functions  (linear  finite  elements) 
h 

v.  corresponding  to  the  triangulation  in  Figure  6.1.  We  set 
h 

K.  ■=  {v.  « V:  v,  > 0 in  D and  v^  = 41  on  3D.  ) . (6.5) 

h h h h — h h 

Since  the  functions  vt  are  linear,  v,  is  non-negative  in  D iff  v is  non-negative 
h h h 

on  D.  u 3Dt  . Since  4>  > 0 on  3d, 
h h — 

K = { v,  tV  : v.  > 0 on  D,  and  v,=  $ on  3D.)  . (6.6) 

h h h h—  h h h 

The  approximation  w^  is  readily  computed  as  is  shown  in  section  7.  Here,  we  de- 
rive an  error  estimate  for  llw-w^H  by  combining  the  ideas  of  Brezzi  and  Sacchi  (to 

appear]  and  Brezzi,  Hager,  and  Raviart  [to  appear] . 


Theorem  6 . 2 


The  piecewise  linear  approximate  solution  w^  exists  and  is  unique.  Furthermore, 


w-whl|  12  = 0(h) 


(6.7) 


Proof:  The  existence  and  uniqueness  of  w.  is  an  immediate  consequence  of  Theorem  6.1 

'*  n 

together  with  the  fact  that  a is  a coercive  bilinear  form. 

2 

We  now  introduce  some  notation.  For  any  two  functions  g^.g,^  L (°)  we  set 


(gl,92)  = J gig2  dxdy 


(6.8) 


2 1 

From  (5.4),  w e H (D)  , so  that  w and  w « II  (D)  and  hence  (see  (4.8)), 

x y 


Lw  = div  exp[f(x)  - g(y)Jgrad  w e L (D) 


(6.9) 


For  any  v„  e H_ (D) 
0 0 


we  thus  have  that 


(6.10) 


l(w,v0)  = / explf(x)  - g(y)]  grad  w grad  vQ  dxdy  , 


= -J  Vq  Lw  dxdy 


= (-Lw,  vQ) 


Finally,  we  note  that 

j(v)  = (e,v)  (6.11) 

where 

e(x,y)  = explf(x)]  . (6.12) 

If  v € K then  v - w e H1 (D)  so  that  using  (6.10)  and  (6.11)  the  variational 

0 

inequality  for  w may  be  written  in  the  equivalent  form 

(-Lw  + e,  v-w)  >0  , (6.13) 

for  all  v e K . 

In  (6.13)  we  may  set  v = w + vQ  for  any  non-negative  vQ  c H^(D),  from  which  we 
conclude  that 

-Lw  + e _>  0 a.e.  in  D . (6.14) 

For  any  E > 0 let  ^ be  a smooth  non-negative  function  which  is  equal  to  1 
at  points  in  D which  are  at  least  a distance  2c  from  3d  and  equal  to  0 at  points 
less  than  a distance  E from  D . Then,  v^.  = V £ w £ H^fD)  • Setting  v = w + v^ 

in  (6.13)  and  letting  E ■+■  0 we  obtain 

(-Lw  + c,  w)  ^ 0 

Setting  v - w - v in  (6.13)  and  letting  E -*■  0 we  obtain 

(-Lw  + e,  -w)  >_  0 . 

Hence, 

(-Lw  + e,  w)  = 0 . (6.15) 

Finally,  by  assumption, 

w > 0 a.e.  in  D . (6.16) 
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W 


Inequalities  (6.14)  and  (6.16)  together  with  equality  (6.15)  constitute  a complemen- 
tarity problem  for  w . 

Since  -a(w.  ,v  -w.)  < (e,  v -w.  ) for  all  v,  e K,  , we  see  that 
hhh—  hh  h h 


a(w-w,  , w-w,  ) = a (w-w.,  w-v,)  + a (w-w,,  v -w,)  , 

h h h h hhh 


a (w-w  , w-v  ) - a (w  , v -w  ) 4 a(w,  v -w  ) , 

h h hhh  h h 


< a(w-w  , w-v  ) 4 (e,  v.-w,)  4 (-Lw,  v.-w.)  , 

— h h h h h h 


= a(w-w.  w-v.)  4 (-Lw  4 e,  v.-w,)  , 

h h h h 


< a2  ||w-wh||1(2  ||w-vh  ||1(2  + (-Lw  4 e,  vh-wh)  . (6.17) 

where  a2  is  the  constant  introduced  in  (5.1).  Using  (6.14)  and  (6.15) , we 
conclude  that 

(-Lw  4 e,  v,  -w.)  = (-Lw  4 e,  v -w)  - (-Lw  4 e,  w,  ) 4 (-Lw  4 e,  w)  , 
h h h h 

= (-Lw4e,  v.  - w)  - (-Lw  4 e,w.)  , 

h h 


< (-Lw  4 e,  v,  - w)  , 
n 

i II  "Lw  + e ||  0f2ll  vh  - w||  0j2 


(6.18) 


2 I 

For  any  v c H (D)  let  v denote  the  piecewise  linear  interpolate  to  v . It 
follows  from  the  work  of  Ciarlet  and  Raviart  [1972]  that  there  is  a constant  C in- 


dependent of  v and  h such  that 


lV  ' vI  Hm,2  - C II  2, 2 h2m'  f°r  «>  = 0.1 


(6.19) 


Next,  we  note  that 


II  W-Wh  II  i ( 2 1 Hw-wI  ll1(2  + 1 1 wI-wh  1 1 1 , 2 * 

2 2 2 

Squaring,  and  using  the  inequality  (a4b)  £ 2(a  4b  ) , we  obtain 

l|w-wh  ||12)2  < ( llw-w1  ||1(2)2  4 ( ||wT-wh  ||12)2 


(6.20) 


Finally,  we  observe  that  w*  - w^  e H^fD)  so  that,  from  the  coercivity  of  a on 
H*(D)  (see  5.3), 


a1<  ||  wI-wh  ||  ^ 2*  2 — a<wI_wh'  wX-wh) 


(6.21) 


Using  (6.19) 


We  can  now  begin  the  final  computations.  We  set  E = ||w-wh  ||  x 2 
(6.20),  and  (6.21)  we  obtain 

• E2  < eye  1 1 w ll2(2>2  h2  + a^-Wj,,  wJ-wh)  , 

= a(wI-w.  # w*-w.)  + 0(h  ) , 

n n 

= a(w*-w,  w -w)  + 2a (w  -w,  w-w^)  + 


Using  (6.19), 


+ a(w-w.  , w-w  ) + 0(h  ) , 

n n 


< a2(  ||wI-w  Il1(2)2  + 2a2  llwT-w  lllj2  II w_wh  II  1,2 


-7  E2  < a,(c||w  II,  y h*  + 2a  C ||w  II  h E + 


+ a (w-w.  , w-w  ) + 0 (h  ) 
h h 


2 . 2 


- 2 


'2,2 


2,2 


+ a (w-w,  , w-w  ) + 0(h  ) , 
h h 


L 


= a (w-w  , w-w  ) + 2C  h E + 0(h  ) , 

h h l 


where  Cj  = C»2  C ||w  ||2^2 


Using  (6.17)  and  (6.18)  with  v - w , 


— E*  < «2(C  ||w||2f2)hB  +J-LW  + e,  w - wh)  + 

+ 2CX  h E + 0(h2) 

I 2 

< 3CX  h E + ll-Lw  + e ||0(2  || w -wh||0(2  + 0(h  ) . 

Using  (6.19)  with  m = 0 , 

-J-  E2  < 3CX  h E + ||  -Lw+e  II 0, 2 C llw  H 2,2  ^ + °<h2>  ' 


3CX  h E + 0(h  ) 
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Thus,  multiplying  through  by  we  have  that 

1 

E2  < 2C2  h E + 0(h2) 

with  C2  = GC^/a^  • Hence, 

(E  - Cjh) 2 < (C2h)2  + 0(h2)  = 0(h2)  , 

so  that 

E - C2h  = 0(h)  , 

and  finally, 

E = 'lW"Wh  Hi, 2 = °(h)  • ° 

Remarks 

1.  There  are  a number  of  interesting  questions  about  the  convergence  of  w.  to 

h 

w: 

(a)  How  fast  does  the  approximate  free  boundary  converge  to  the  true  free 
boundary?  In  this  connection  see  Brezzi,  Hager,  and  Raviart  [to  appear, 
p.  21]  and  Brezzi  and  Sacchi  [to  appear,  p.  9]. 

CO 

(b)  Can  one  obtain  an  L estimate  for  the  error?  In  this  connection  see 
Baiocchi  [1976a]. 

2.  We  draw  attention  to  a number  of  related  references  on  the  numerical  solution 
of  variational  inequalities:  Falk  [1974];  Glowinski  [1976];  Glowinski , Lions,  and 
Tremolieres  [1976];  Mosco  and  Strang  [1974];  Hager  [1976];  Hager  and  Strang  [1975]. 

3.  The  functions  v.  must  satisfy,  at  least  approximately,  the  boundary  conditions 

h 

v = $ on  3d  and  the  inequality  restraints  v ^ 0 a.e.  in  D . By  formulating  the 
problem  in  terms  of  v *=  v - 4>,  the  boundary  conditions  take  the  simple  form  v = 0 on 
3d,  but  the  inequality  constraints  take  the  more  complicated  form  v >_  -4>  . It  was, 

therefore,  decided  to  retain  the  formulation  in  terms  of  v . 
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7.  A numerical  example 


As  an  example  we  consider  the  specific  geometry  shown  in  Figure  7.1,  which  was 
chosen  because  it  had  previously  been  considered  by  several  authors. 


Fiqure  7.1:  A numerical  example  (r  = 4.8,  R = 76.8,  h = 12,  H = 48, 



m = 8,  n = 12 . ) 

Because  the  solution  changes  most,  rapidly  near  the  well,  the  subdivisions  were 
taken  to  be  uniform  in  the  y-dircction  and  logarithmic  in  the  x-direction.  If  n and 
m denote  the  number  of  subdivisions  in  the  x-  and  y-directions,  the  coordinates  of  the 
gridpoints  were  given  by 

Vj  = j H/m,  0 £ j £ m , 

= r expKi/i’J  ln(R/r)],  0 <_  i <_  n 

The  integer  m was  always  chosen  to  be  a multiple  of  4 so  that  the  corner  D was  a 
gridpoint;  this  was  advisable  since  w is  not  smooth  at  the  corner  D . 
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The  permeability  K.  was  taken  to  be  one  so  that  f (x)  = In 


and  g(y)  = 0 


From  (4.23),  (4.24),  and  (6.2), 


J(v)  = / x[v2  + v2  + 2vjdxdy 
D * y 


(7.1) 


From  (4.12),  (4.13),  (4.14),  (4.15),  and  (4.18),  the  boundary  values  $ are 

4>(x,H)  =0,  on  AF  , 

<t>(R,y)  = (H-y)  2/2,  on  AB  , 


>(r,y)  = (h  -y)  /2,  on  CD  , 
w 


(7.2) 


4>(r,y)  = 0,  on  DF  , 


h2  ln(R/x)  + II2  ln(x/r) 


4>(x,0) 


2 1 n ( K/r ) 


, on  BC 


For  v c K,  let  V = {V  ) denote  the  N = (m+1)  * (n+1)  vector  such  that 
h h ij 


V. . = v (x. ,y  . ) 
i]  h l j 


(7.3) 


Then 


J(v,_)  = VTAV  + 2bTV 
n 


(7.4) 


where  A - { A ( i , j ; i,j)}  is  an  N X N matrix  and  b = {b(i,j) } is  an  N-vector.  The 
matrix  A and  vector  b are  best  obtained  by  computing 

.2  . . .2 


W = /„  X,(Vx)  + ‘"h.y’  + 2Yh)dxdy 


(7.5) 


= VT  ArV  + 2b^V  , say, 


(7.6) 


for  a grid  rectangle 

R = {(x,y):  x.  < x < x.+1,  y . < y < yj+1)  , (7.7) 

and  then  summing  over  all  grid  rectangles. 

Consider  such  a rectangle  R which  we  divide  into  two  triangles,  L and  U . 
(See  Figure  7.2) . 
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Then 


On  L: 


On  U: 


v.  = (V.  , . - V.  . )/Ax  , 
h,x  i+l  #3  1,3 


v.  - (V.  - V.  , )/Ay  , 

h,y  1,3  + 1 i,3 


v = V.  . + (x-x.)v.  + (y-y.)v 

n 1,3  a h , x 3h,y 


v-  = (V.  - V.  . . )/Ax  , 

h,x  l+l, 3+1  1,3+1 

v — (V  - V . )/Ay  , 

h,y  i+l, j+l  1+1,3 


vh  = vi+i,j+i  4 <x-xi+i)vh,x  4 vh,y 


Let  X and  X be  the  x coordinates  of  the  centroids  of  L and  U : 
L U 

xr  = x.  + Ax/3  , 

L 1 

x^  = x^  + 2Ax/3 


(7.8) 


(7.9) 


(7.10) 


By  direct  computation,  the  non-zero  components  of  A are: 


AR(i,j;i,j)  = - Ax  Ay  xl(Ax  2 + Ay  2) 

1 -2  -2 

AR(i+l, j+lji+1, j+l)=  — Ax  Ay  xy(Ax  + Ay  ) , 

1 -2  1 -2 

A ( i + 1 , j ; i + 1 , j ) = — Ax  Ay  x Ax  + — Ax  Ay  x Ay  , (7.11) 

R 2 L 2 U 

1 -21  -2 
AR(i, j+lii, j+l)  = 2 Ay  XL  Ay  + 2 Ax  Ay  xu  ^ ' 

1 -2 

AR (i,  j ,•  i+1 , j)  = --  Ax  Ay  xL  Ax  , 

1 -2 

AR ( i , j j i , j + 1 1 = --  Ax  Ay  xL  Ay 

Ar( i,j  + l;i  + l,j  + l>=  Ax  Ay  xy  Ax  2 , 

1 -2 

AR(i  + l,  j ; i + 1 , j + l)  = - - Ax  Ay  xy  Ay  , 

A (i+1, j ; i , j ) = A , 

K K 

AR(i,j+l;i,j)  = AR(i.j;i,j+l). 

A (i+1 , j+l ; i , j+l)  = A (i,j+l;i+l, j+l)  , 

K K 

A (i+1, j + l i i+1 , j)  = AR( i+1 , j ; i + 1 , j + l)  • 

The  non-zero  components  of  b are : 

K 

..  1 . . 1 A A ,*L  j.  2Ax  , - 

Vl,j>  = - Ax  Ay  xL  - - Ax  Ay  [y  + — ) 


bR(i+l,j+l) 


The  problem 


1 . . 1 . . , U 2Ax  , 

= 2 Ax  Ay  Xu  - - Ax  Ay  I— - - — ] 


1 A A f U Ax 

- I Ax  av(T  + 3£  1 • 


T T 

Minimize:  V AV  + 2b  V 


Subject  to:  V. . = ♦. . on  3d  , 
13  13  h 


(7.13) 


V.  . > 0 , 

13  “ 

is  a quadratic  programming  problem  which  is  equivalent  to  a complementarity  problem. 

Many  algorithm:',  are  available  for  the  solution  of  this  problem:  see  Cottle  11974,  1974a), 
Cottle,  Golub,  and  Sachcr  [1974],  Cottle  and  Pang  [1976),  Cottle  and  Sacher  [1973]. 

Here,  we  use  a valiant  of  S.O.R.  (systematic  overrolaxat jon)  to  solve  (7.13). 

Given  an  initial  guess  satisfying  the  boundary  conditions,  we  generate  a sequence 


of  approximations 

v(k) 

using  th< 

following  AT.GOT,  segmen 

For 

i :==  1 

step 

i unti 1 

n-1  do 

For 

j :=  1 

step 

1 until 

m-1  do 

vtl  - 

[b  (3 , j)  - 

V(i,j+1)*  A ( i , ;i ; i , j H 1 ) - 

(7.14) 


- V (i-t  1,  j ) * A ( i , j ; i+1 , j ) - 

- V ( i - 1 , j ) * A (i , j ; i-1 , j)  - 

- V(i,j-1)*  A ( i , j ; i , j- 1 ) 1 /A (i , j ; i , j ) ; 

vt2  = V(i,j)  + W*  (vtl  - V (i , j ) ) ; 

v(i, j)  - if  vL2  £ 0 then  0 else  vt2  ; 
end; 

Here,  (0  is  an  overrelaxation  parameter  which  must,  be  chosen  to  optimize  the  rate  of 
convergence . 

The  algorithm  (7.14)  is  known  to  converge.  The  algorithm,  and  related  algorithms 
have  been  considered  by  many  authors:  Hildreth  [1957];  Merzljakov  [1962);  Fridman  and 
Chcrnina  [1967];  Martinet  [1967);  Durand  [1969/1969,  1972);  Glowinski  [1971,1973]; 
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Cryor  [1971);  Comincioli  [1971);  Miellou  [1971,  1971a,  3 972) ; Luong  [1973);  Martinet 
and  Am; lender  [1974);  Eckhardt  (1974);  Mangasarian  [1976]. 

To. choose  the  parameter  W we  proceeded  as  follows.  For  the  case  m = 16,  n = 24 
computations  were  made  with  several  values  of  w and  it  was  found  that  the  optimum 
value  of  u was  approximately  1*7  . For  the  general  case  we  set 

U = 2 

1+sxnm/nf  ict) 

where 

nfict  = *8590  [ (nt 1) (ra+1) ] 1/2  . 

The  expression  for  to  is  a modification  of  the  theoretical  optimim  value  for  S.O.R. 

on  a square  (Varga  [1962,  p.  203)).  The  constant  *8590  was  chosen  so  that  to  = 1*7  when 

m = 16  and  n = 24. 

The  computations  presented  no  difficulties.  The  solution  of  the  smallest  problem 
is  given  in  Table  7.1. 

In  Table  7.1  the  position  of  the  approximate  free  boundary  is  shown  by  the  first 

zero  term  in  each  column.  The  approximate  solution  is  identically  zero  on  the  vertical 

line  x = r so  that  it  is  not  possible  to  determine  the  height  h^  at  which  the  free 

boundary  intersects  the  well.  As  an  approximation  to  h we  take  the  height  h’  of 

s s 

the  free  boundary  at  the  vertical  gridline  adjacent  to  the  well.  For  example,  from 

Table  7.1  we  obtain  h'  = 36. 

s 

In  Table  7.2  the  values  of  h^  for  a sequence  of  decreasing  grid  lengths  are  given. 
The  iterations  were  terminated  when, 

||vm  _ v10(a-l)||  <1q-6  _ 

00 

In  judging  this  accuracy  criterion,  it  should  be  remembered  that  ||v  ||  = 1152,  so  that 

-9 

the  relative  error  is  10 
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r w 


i 


ro 

n 

h* 

s 

number  of 
iterations 

4 

6 

36.00 

20 

8 

12 

30.00 

40 

16 

24 

30.00 

70 

32 

48 

30.00 

120 

64 

96 

30.00 

230 

Table  7.2:  Values  of  h' 


The  method  of  determining  u>  seemed  satisfactory,  although  the  ratios 

||Vi°k  - v10(k-% 

||v1°(k-l)_v10(k-2)||2 


oscillated  so  as  to  suggest  that  the  dominant  eigenvalue  of  the  iteration  was  complex  and 
hence  that  (from  the  theory  of  S.O.R.)  reducing  (U  would  improve  the  convergence. 

For  comparison,  we  compare  in  Table  7.3  the  values  of  hg  obtained  by  different 
authors.  With  the  exception  of  the  present  computation,  all  the  results  are  presented 
graphically  so  that  we  have  had  to  estimate  h from  graphs. 


Author 

Method 

h 

s 

Hall  11955,  p.  29] 

trial- free-boundary; 
finite  differences 

34.0 

Taylor  and  Luthin  [1969] 

time-dependent;  finite 
differences 

34.0 

Neuman  and  Witherspoon 
[1970] 

trial  free-boundary; 
finite  elements 

30.0 

Neuman  and  Witherspoon 
[1971,  p.  620] 

time-dependent;  finite 
differences 

30.0 

Present 

Variational  inequalities 

30.0 

Table  7.3:  Computed  values  of  hs 
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The  differences  in  Table  7.3  may  be  explained  by  the  fact  that  the  physical  assumptions 
differed:  Hall  [1955)  assumed  capillarity  and  a lined  well;  Taylor  and  Luthin  [1969] 

assumed  partially  saturated  flow;  and  Neuman  and  Witherspoon  [1970,  1971)  made  the  same 
assumptions  as  in  the  present  paper. 

Finally,  in  Table  7.4  we  give  the  computed  values  of  w at  a typical  point  x = 19.2 
y = 24  for  different  values  of  m and  n . We  also  give  the  differences  between  succes- 
sive approximations  and  the  ratios  of  successive  differences.  It  can  be  seen  that  the 
results  are  consistent  with  the  hypothesis  that 

2 

w-w.  = 0(h  ) . 
h 

In  Theorem  6.2  we,  of  course,  only  proved  that 


= 0(h)  . 


m 

n 

Wh 

Aw. 

n 

ratio 

4 

6 

81.333841 

0 

12 

83.478767 

2.144926 

16 

24 

84.150541 

.671774 

3.193 

32 

48 

84.298332 

.147791 

4.545 

64 

96 

84.337269 

.038930 

3.796 

Table 

7.4:  Values  of  w, 
h 

at  x = 19.2, 

y = 24.0 

In  conclusion  we  draw  attention  to  some  related  work: 

Numerical  solution  of  variational  inequalities  for  porous  flow 
Comincioli  [1974,  1974a,  1974b,  1975], 

Comincioli,  Guerri , and  Volpi  [1971], 

Cottle  [1974]. 

Numerical  sol ution  of  axisymmetr;  <-  wel  1 problems 
Babbit  and  Caldwell  [1948]; 

Yang  [1949];  Boulton  [1951];  Kashef,  Touloukliam,  and  Fadum  [1952); 

Kashef  [1953];  Schmidt  [1956];  Murray  [I960];  Kirkham  [1964];  Taylor  [1966]; 
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Taylor  and  Brown  [1967]  (see  also  Kealy  and  Busch  (1971);  Herbert  [1968]; 
Mauersberger  [1967,  1968,  1968a,  1968b,  1968c];  France,  Parekh,  Peters,  and 
Taylor  [1971]. 


Appendix  A : An  alternative  numerical  approach 

Brezzi  and  Sacchi  [to  appear]  have  givon  a different  error  analysis  for  the  case  of 
plane  flow  through  a rectangular  dam. 

The  approach  of  Brezzi  and  Sacchi  requires  that  be  chosen  so  that  for  every 

v.  c K,  there  is  a v c K with  v < v,  . In  this  appendix  we  show  how  one  can  choose 
n n — h 

so  that  this  condition  is  satisfied. 

In  section  6 the  approximation  was  obtained  by  triangulating  D and  approximating 
w by  piecewise  linear  functions,  the  boundary  conditions  b^ing  satisfied  approximately 
by  interpolation  at  the  nodes.  However,  in  the  special  case  when 

f (x) 

e = x , 

we  have  from  (7.2)  that,  on  BC  , 


so  that 


$(x,0)  = [t(R,0)  ln(x/r)  + <i’(r,0)  ln(  R/x ) ] /In (R/r ) , 


4>  (x,  0)  = -[<>(R,0)  - ‘Mr,0]/x  In  (R/r)  , 


< 0 . 


Thus,  <!>  is  concave  on  BC  and  any  linear  interpolate  v lies  below  any  v c K . 
fore,  even  in  this  case,  the  condition  of  Brezzi  and  Sacchi  is  not  satisfied. 

To  overcome  this  difficulty,  we  introduce  now  coordinates  a and  fl  defined  by 

y 


There- 


a = / e 
r 


-f  (t) 


dt,  P = / eg^  dt  . 


(A.  1) 


The  rectangle  r>  is  transformed  into  the  rectangle  D’ : 

R 


0 < a < / 


-f  (t) 


dt  - aR,  say, 


(A. 2) 


0 < P < / eg*"  dt  = B , say. 


We  have  that 


w = w a 
x a 


= vj  e 
* a 


-f(x) 


w = w„  3 = wa  e 

y p y 3 


g(y) 


(a. 3) 
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The  Jacobian  of  the  transformation  is 


9x 

9x 

9a 

96 

h. 

h. 

9a 

96 

9a 

9a 

9x 

9y 

96 

96 

9x 

9y 

-f  (X) 

0 

0 

e9(y) 

-1 


-1 


f(x)  - g (y ) 


(A. 4) 


Thus  the  minimization  problem  (6.2)  is  equivalent  to  the  problem: 

min  J' (v' ) , 

v‘cK* 


where 


(A. 5) 


K'  = {v*  c H1(D')  : v'  - 4>’  £ H*(D') 

and  v'  > 0 a.e.  in  D')  > 


(A. 6) 


J'  (v*)  = a'  ( v 1 ,v’)  +2  j'  (v’)  , 


(A. 7) 


(w'jV1)  = / / (e  w'  v'  + e2f*X^  wi  vi)dad6  , 

D’  a a ts  P 


(4.8) 


jMV)  = / / e2f(x>-9<y>  vi  dadB  . 

D' 


(A. 9) 


We  note  that  f(x)  and  g(y)  are  continuous  so  that  the  mapping  (x,y)  (a,B)  is 

l'smooth  and  (D)  is  mapped  homeomorphically  onto  H1(D’)  (Adams  [1975,  p.  63))- 


-37- 


We  can  obtain  an  appi-oximate  solution  wR  by  proceeding  as  in  section  6.  The 

rectangle  D'  is  triangulated  in  the  same  way  as  shown  in  Figure  6.1.  The  subdivisions 

are  not  uniform,  but  the  ratio  (interval  length)/h  is  bounded  both  above  and  below.  The 

interior  gridpoints  are  denoted  by  D'  and  the  boundary  gridpoints  by  3d'  . The  space 

h h 

of  piecewise  linear  functions  on  D'  is  denoted  by  V'  . 

h h 

The  approximation  w'  is  then  given  by: 
h 


J'(w')  = Min  J’(v')  , 

v;  c k;  h 
h h 


(A. 10) 


Kj  = {v'  c V'  : v'  > 0 in  D' 
h h h h _ h 


and  v'  = C1'  on  ?D'} 
h h 


Lemma  A . 1 


If  v^  e ‘J’1  - v^  < 0 on  3d'  . 


Proof:  On  F^  we  have  from  (4.16)  that 


(fCo^O)  = — fct  <I>(r,o)  + (a  -a)  <f  ( i , 0 ) J , 
aR  R 


(A. 11) 


so  that  I'1  is  linear  and  hence  v*  - 4>'  on  1''  . 

h a 

On  we  have  from  (4.13)  that 

y 

4>'  (a  ,B)  = Constant  - / eg^(H-t)dt  , 

R 0 

^ y 

= Constant  - Up  + / e9  ^ t dt  . 

0 

Let  (a  ,6  ) and  (a  ,6  ) be  two  adjacent  nodes  on  F ' corresponding  to  (R,y  ) and 
R j R 2 i-  1 

(R,y2)  on  Fx  . We  may  assume  that  B.,  > B^  . Then,  for  B^  £ B _£  B2> 


!>’(«„, B)  - v.  (a  ,B)  = / e9(t)(t  - c)dt 

R h K ' 

y^ 


= E(y)  , say  , 


where  the  constant  c is  such  that  E(y2>  = 0 • Wo  conclude  that  y^  < c < y2 
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Thus,  E(y)  decreases  monotonically  for  y £ y <_  c and  increases  monotonically  for 
c 1 V £ Y2>  so  that  E(y)  £ 0 . 

Similar  arguments  show  that  4>'  - < 0 on  ]’£,  while  <J>'  - = 0 on  Tj  u Tj.  - D 

Theorem  A. 2 

The  approximate  solution  w^  exists  and  is  unique.  Furthermore, 


1,2 


= 0(h)  . 


Proof : The  existence  and  uniqueness  of  w^  is  a consequence  of  the  fact  that  the  bi- 

linear form  a'  is  coercive. 

Remembering  that 


-a* (w' ,v'-w' ) £ 

j* (v'-w'J, 

, for 

all 

v'  e K'  , 

and 

-a‘ (w‘ ,v'-w')  < 
h h h — 

jMv'-w') 

, for 

all 

Vhe  K ' 

we  see  that  for  all  v'  c K*  and 

Vh  f "h  ' 

a* (w'-w’ , w'-w/) 
h h 

= -a'  (w' 

,„'-W 

) - 

a' (Wh'W'-Wh>  ' 

= -la' (w* ,v'-w')  + a' (w* ,w'-v' ) ] - 

h 

-la'  (w ' , v ' - w ' ) + a'(w',w'-v')J  , 
h h ■ h h h 

< j'(v'-w')  - a' (w', w'-v')  + 

— h 


+ j'  (v'-wM  - a'  (w'  ,w'-v')  . 

h h h h 

= -j' (w'-v'+w'-v' ) - a' (w’, w'-v')  - 
n n n 

-a1 (w * , w'-v')  . 
n n 


Adding  and  subtracting  the  term  a' (w',w'-v')  from  the  right  hand  side  we  obtain 

h 

(following  Brezzi  and  Sacchi (to  appear]), 

a' (w'-w',  w'-w')  < -j' (w'-v'  + w'-v')  - 
h h — h h 


-a'(w'#  w'-v'  + w ' - v')  + 
n n 


r 


■ 


a'fW-w'.w'-w')  < -a'(w',»'-v')]  + 


+ [-j'(w’-v')  - a 1 (w’ , w'-v' ) J + 
" h 

+ a ' (w1 -w' , w' -v' ) , 

h h 


(A. 12) 


for  all  v^  c and  v'  e k'  . 

The  remainder  of  the  proof  is  very  similar  to  the  proof  of  Theorem  G.2. 

2 

By  (5.4),  w c H (D)  . We  have  assumed  that  f and  g arc  continuously  differ- 
entiable so  that  the  mapping  (x,y)  - <a,P)  is  2-smooth  and  1!2<D)  is  mapped  homeo- 
2 

morphically  onto  H (D * ) (Adams  (1975,  p.  63)).  Thus, 


w*  c 11  (O')  . 

Let  w'  denote  the  piecewise  linear  interpolant  to  w'  . From  Theorem  5 of 
Ciarlet  and  Ravi  art  11972]  wo  conclude  that 

IJw  - W*1  11^,2  < c J|w  \\2'?  h2_I",  for  m = 0, 1 . 

Using  the  came  arguments  and  notation  as  in  Theorem  6.2  we  obtain  tliat 

a' (w* ,v^)  = - (L'w1 ,v^)  , 

for  all  v^  c 11*  (D*),  where 

L'w'  = h <e'29(y>  wa>  + -k  (e2f<X)  w6>  £ l2(d,)  • 

Let  j‘  (V)'  = (e*  ,V)  , 


(A. 13) 


(A. 14) 


(A. 15) 


where  e' (a,B)  = exp[2f(x)  - g(y))  . 

Then  the  identities  (A. 35)  and  (A. 37)  can  be  used  to  manipulate  (A. 12).  We  obtain 

a’  (w’-w^.w’-WjM  <_  [— j 1 (w1  “V* ) - a' (w>  ,w'-v|  )]  + 

+ l-j'W-v')  - a'  (w'  ,w'-v' ) ) + 

h h 

+ a' (w1 -w' , w'-v' ) , 

h h 
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(A.16) 
(A. 17) 
(A.1G) 


Ij'(w'-V^)  + j'(w^-v’)) 


[a* (w' ,w'-v' ) + a'(w',  w'-v')]  + 

n n 


♦ a* (w'-w' , w'-v* ) , 
h h 


= - I (e'  / W -v* ) + (e',w'-v')l  - 
n n 


■ [ (-Lw' ,w' -v' ) + (-Lw' ,w'-v' ) ] + 

n n 


♦ a' (w' -w' , w'-v' ) , 
n h 


*=  -(-Lw'  + e' , w’-v' ) - 
n 

-(-Lw'+e'#  w'-v')  + 
n 

+ a ' (w'-w' , w'-v')  (A. 19) 

h n 

where  in  the  last  step  but  one  we  have  used  the  fact  that  w'-v'  and  w'-v*  belonq  to  H3 (D' ) . 

h h ^ 0 

In  (A.  19)  we  now  set  v'  = w'  . By  Lemma  A.l  we  have  that  4>'  - v'  < 0 on  3d' 

h h — 

so  that  v’  - <_  0 on  3d'  for  all  v'  e K'  . Since  is  non-negative , we  may 

choose  v'  e K'  so  that  v'  - v^  0 on  D'  . Analogously  to  (6.14)  wc  can  show  that 

-L'w'  + e'  ^ 0 a.c.  in  D'  . 

Thus,  in  (A. 19), 

-(-L'w'+e',  w'-v')  = -(-Lw'+e',  w,:l-v')  < 0 . 

n — 


Noting  (A. 13)  and  (A. 14)  we  thus  conclude  that 


a' (w'-w^,  w'-w^)  < ll-Lw'+e'  ||1<0  C |(w’  ||2  2 h + 

+ tt2  HW'-Wh,ll,2  C l'W'  H 2,2  h ' 


ci  ,lw''wh  Hi, 2 h + 0(h  > • 


(A. 20) 


where  a 2 is  the  coercivity  constant  for  a'  and 


c;  = liw' il2,2 c • 


As  in  (6.20)  and  (6.21)  wo  have  that 
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Appendix  P : The  computer  program 

The  computer  program  .used  to  obtain  the  numerical  results  quoted  in  Section-7  is  listed 
below.  Two  minor  remarks  are  perhaps  necessary: 

1.  The  main  program  is  written  in  ALGOL.  Since  the  ALGOL  compiler  available  to  us 
does  not  optimize,  the  inner  S.O. R.  loop  is  executed  by  a FORTRAN  subroutine. 

2.  The  computations  were  performed  in  double  precision  so  that  the  asymptotic  behavior 


of  the  error,  as  shown  in  Table  7.4,  would  not  be  contaminated  by  round-off  errors. 


•NUALG.ISZX 

BEGIN 

COMMENT  **** ************ **•*«•**•**»•**•*»*« ***************  ********* 

PROGRAM  FOR  THE  numerical  SOLUTION  , BY  MEANS  OF  THE 

variational  method  proposed  by  baiocchi  et  al,, 

OF  THE  FREE  BOUNDARY  PROBLEM  RELATED  TO  THE  STATIONARY  FLOW 
THROUGH  A FULLY  PENETRATING  WELL 


' 


VARIABLE  names  used  in  the  program  INCLUDE  the  FOLLOWING  l 


«Y1*  IS  THE  height  OF  THE  WATER  AT  X a B 

*Y 2*  IS  THE  HEIGHT  OF  THE  WATER  AT  X a A 

«B*  IS  THE  RADIUS  OF  THC  CATCHMENT  AREA  OF  THE  WELL 

tA*  IS  THE  RADIUS  OF  THE  wfll 

AND  *N*  APE  THE  NUMBERS  OF  POINTS  OF  SUBDIVISION  OF  THE 
6IDES  OF  THE  RECTANGLE  *R», 

IOMEGA*  IS  THE  RELAXATION  PARAMETER, 

*EPS*  IS  USED  IN  THE  TEST  FOR  STOPPING  THE  ITERATION, 

ARRAY  *U*  SHALL  CONTAIN  THE  FUNCTION  VALUES  OF  THE  DISCRETE 
SOLUTION. 

EXTERNAL  FORTRAN  PROCEDURE  ITERATf 
COMMON  COEFO  ( R E A L 2 ARRAY  CO  (0 1 4« , 0 1 32 > » ) » 

COMMON  COEFJ  (REAL2  ARRAY  Cl  ( 0 t UP , 0 I 32 ) f 1 | 

COMMON  C0EF2  (REAL2  ARRAY  C2  (0  I 4 A, 0 | 32  ) » ) | 

COMMON  COEF3  (REAL2  ARKAY  C3  ( 0 J 4B  , 0 I 32 ) > ) | 

COMMON  UNK  (REAL2  ARRAY  U ( 0 I 48 , 0 I 3? ) » ) » 

COMMON  PAR  (INTEGER  I m a X , N, H, MOD  I T j RE AL2  OMEGA, TEST, TEST2 >) I 
REAL2  ARRAY  R (0  » 48)> 

REAL2  ARRAY  S (0  I 32)> 

REAL2  DELTAX,DELTAY» 

REAL2  NFICTr 

Rf AL2  A, VI, Y?,EPS, HI, H2i 
REAL2  A1,B1,A1B1| 

REAL2  B,AB,DENOM» 

REAL2  MU, LAMBDA, PJ» 

RE AL2  TE8T2P,RATI0j 
RCAL2  T!,T2,R2,RS» 

REAL2  T1 P , P2P , RSP > 

INTEGER  !,J,K,ITER,LB,UB> 

INTEGER  NOCOLS,NP> 

COMMENT  a***************************** ***********fF**e*** *«***»•*«* 
FORMAT  FMTl  (El, X10, 'NUMBER  OF  ITERATIONS  a f,I<t,A1.0, 
Xn,7(R15,8,Xn,A1.0, 

Xll,7(15('-'),xn,A1.0), 

FMT2  (X2,06t4,X3,7(R15,6,Xn,Al,0), 

FMTU  (X5, ' I TER  « ',!3,'  TEST  a ',R1U,7,»  L2  NORM  p »,R14t7, 
' RATIO  a ',Rlfl,7,AJ,0), 

FMTO  (E 1 , X5 , » RW  a ',R12,U,*  RE  a ' , R1  ?.  4 , A 1 , 0 , 

X5, ' HW  a ',R12,a,i  HE  a ' , R12.0, Al ,0, 

X5 , ' NX  a ',  112  , • NY  a H2  ,A1,0, 

X5, ' EPS  a ' ,R12,4,  ' OM  a ' , Rl 2 , 4 , A 1 , 0 ) , 
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FMT5  (X5, '***', X2, 'TEST  a ',  Rla,7,  '***',*1.0)1 

COMMENT  ...........  GIVEN  (HUNT JTIES  , 

READ  (CARDS,M,N,A,AB,Yl,Y2,fP3)| 

I m 4 X |a  46  * If 

P!  is  a.ORRO  * ARCTAN  (1#0440)| 

*1  l«  LN  ( 4 / A ) | 

B l«  A ♦ AP | 

B1  |s  UN  (B/A) | 

A 1 B 1 | a B l • A 1 > 

HI  | a AlBl/Nj 
HZ  | a Y1  / M, 

COMMENT  *H l * ANO  *H2*  ARE  THE  STEPS  1 7ES  FOR  THE  OISCRETIZATIONf 
NFICT  | a 0.65B0440  * SORT  (1,0440  * (N*l)  * (M*l))| 

OMEGA  |s  2,0440/(1  ,0440  ♦ SIN  (PI/NFICT))» 

COMMENT  COMPUTE  OVERRELAXATION  PARAMETER  OMEGA  » 

FOR  I |a  0 STEP  1 UNTIL  N DO 
R ( I ) 1=  A *F.XP  ( I*H1  ) j 
FOR  J |=  0 STEP  1 UNTIL  M 00 
S (J)  l a J*Y1/M| 

WRITE  (PRINTER, FM TO, A ,B  , Y2 , Y1 , N, M, EPS , 0MCG A) > 

ITER  is  0 1 

test  t = i,om» 

TEST2  1=  0,04*0 1 
TEST2P  is  1 , 0440 | 

COMMENT  **************************************************************** 
INITIALIZE  *U(I, J)* 

ON  THE  BOUNDARY  OF  THE  BECTANGLE  *R*  , 11(1, J)  IS  GIVEN  BY 
THE  FUNCTION  G(X,Y)  (CF,  BAICCCHI  ET  AL,) 

EVERYWHERE  ELSE  WE  SIMPLY  SET  IT  TO  SOMf.  CONSTANT  GE  0 j 

COMMENT  ———.—ON  THE  SEGMENT  (A,B)  — — > 

FOR  I 1=  0 STEP  1 UNTIL  N DO 

U(I,0)  is  0,54S0*CY1**2*IN(P(I)/A)+Y2**2*LNCB/R(I)))/H1» 

COMMENT  .—.——-on  the  segment  (6,0  ----------- | 

FOR  J is  l STEP  1 UNTIL  H DO 

U(N,J)  |b  0,5440  * (Y1  - Y1  * J / H)  **2 | 

COMMENT  ———ON  THE  SEGMENT  ( A , D ) .....  ..... »t 

K fa  ENTIER  (Y2  / H2)  | 

WRITE  (PRINTER,  <<I3  ,Al»,K)j 

FOR  J | a 1 STEP  l UNTIL  K DO 

U <0,J)  is  0,5440  * ( Y2  - Yi  * J/  M)  **2  | 

COMMENT  — INITIALIZE  U (I,J)  EVERYWHERE  BY  LINEAR  INTERPOLATION  — — | 
FOR  I I*  1 STEP  1 UNTIL  N-J  DO 

FOR  J I a 1 STEP  1 UNTIL  M-l  DO 

BEGIN 

LAMBDA  |«  (R(I)  • R(0))  / ( R ( N ) . R ( 0 ) ) | 

U (I,J)  la  U(N,J)  « LAMBDA  ♦ U(0,J)  * (1,0440  • LAM8DA)| 

END) 

Comment  ************************************************************** | 
FOR  I is  0 STEP  1 UNTIL  N . 1 OP 

FOR  J |iO  STEP  1 UNTIL  M - I DO 

BEGIN 

OELTAX  I a R ( I + 1 ) . R(I)| 

DELTAY  la  S(J*i)  - 8 ( J) | 

COMMENT  (LOWER  TRIANGLE)  — — — — — — — — — — 


Til*  (R(I)  ♦ 1 * PELTAX  / 3.0*00)  * DELTAX  * PELTAY  / 2, 
T2  »■  DELTAX  **?  * DELTAY  /36.0UIP, 

R2  |a  PELT  AX  * T2  ♦ (R(J)  ♦ l • DELTAX  / 3, 0**0)  * TI, 


RS  ,a 

m 

PELTAY  * T?  / 

2 ♦ 

(S(J)  ♦ 

1 * PELTAY  / 3 , 0 f-  X 0 ) * TI, 

C0(I 

) 

!=COd 

,J  ) 

♦ T 1 / 

PEI.  T A X A*?, 

Old 

) 

1=01(1 

,J  ) 

• T 1 / 

peltax  **?, 

00(1 

> J 

) 

taeod 

,J  ) 

♦ T 1 / 

PFLTAY  * *2 , 

C2  (I 

) 

l=CS(I 

, J ) 

- T 1 / 

PELTAY  **?, 

C 3 c I 

) 

t =03 ( I 

,J  ) 

♦ Ti 

- (RS  - R(I)*Tl)/  PELTAX 

-(PS 

- S(J)*T1)/  PELTAY, 

COd+l 

) 

1=00(1*1 

r J ) 

♦ T 1 / 

DELTA):  **? , 

03  d ♦ 1 

, J 

) 

1 aO  3 ( I ♦ 1 

, J ) 

♦ ( R2 

- R(I)*Ti)/  PELTAX, 

COd 

,J+n 

» *C0 ( I 

, J+l) 

♦ TI/ 

PELTAY  **?, 

C 3 ( I 

, J+n 

1 =C 3 ( I 

, J+l  ) 

♦ ( RS 

- S(J)aT1)/  deltay, 

comment 

(UPPER 

TRIANGLE) 

m 

TiPra  (R(I)  +2  * PELT  AX  / 3.0SCP)  * DELTA*  t PELTAY  / ?, 

RPPja  PELTAY  * T2  + (R(I)  <2  * DFLTAX  / 3,CIU0)  * TIP, 

PSP  | a - DELTAY  * T2  / ? ♦ (S(J)  +2  * PEl.TAY  / 3.0P-S.0)  * TIP; 

C 0 (I ♦ 1 , J « l ) tsCOf  HI, .1*1)  ♦TIP/  PEL  T A X **?, 

ItCOdd  , J*1 ) + T1P/  PELTAY  **?, 

I =C3  (!*I|.IM)  ♦ T 1 P+  (R?p-  P (1  + 1)  * T1 P) /DtLT  AX 

♦ (PSP-  S ( J*1  ) *T  1 1’)  / PE  I T A Y , 


co(i4j,  j+n 

C3(I  + 1 » J+l ) 


TIP 

TIP 


co  (i  ,j+n 
Cl  ( I ,J*i) 
c 3 c i ,jm) 

CO  ( 1 + 1 , J ) 

C 2 C T ♦ l , J ) 

C 3 ( I ■*  J ft!  ) 

F wo  r 

COH^CNT  *♦*-♦*•***««*«**********•*«  tA**************************  A.  **«•**! 

FOP  ITEP  i=  ITER  + 1 WHILE  TEST  CEP  EPS  AN|>  1 T(  R LEO  <>00  PI) 

BCCI‘J 


I =CO  C I #.1  + 1) 
x * c i ci  ,j+n 
t sC3  ( i ,Jtn 
I s C O ( I ♦ 1 » >1  ) 

J sC2  d + 1 , J ) 
I =C3 (I ♦! , J ) 


/ PELTAX*r.2j 
/ PEI  T A X * * 2 j 
-(P2P  - R(H15*TlF) 
♦ tip  / DELTA V»«2, 

- TIP  / PEL  TA  Y*<>2  I 
-(RSP  - SCJ+1 ) *T 1 p ) 


/ DELTAX, 


/ PELTAY, 


MPDIT  :p  “OD  (ITER,  10), 

ITERAT  , 

IE  HOPIT  EQL  0 THCN 
BEGIN 

TEST?  | e SORT  (TEST2/C (M-l)* (v-i ) ) ) , 

RATIO  la  (TE. ST?/TFST2f’ !>♦♦©,  1 , 

WRITE  (PRINTER, EWTR, ITER, TE3T, TESTS, RATIO), 

T C S T ? P i o TEST?, 

TESTS  is  0 1 0 it  0 , 

END, 

ENP, 

COMMENT  #****#**#*, :***********  ti  **************  ft ************ 
UB  I = 0, 

LB  1 = 1, 

ITCR  fa  ITTR  - 1 , 

N0C0L2  tc  N ♦ 8, 

FOR  NOCPLS  ia  NOCPLS  - 7 WHILE  NOCPLS  GTR  0 PD 
BEGIN 

HP  ,a  MIN  ( NOCOL 9,7), 

UB  t = 1>B  ♦ NP, 

WRITE  (PRINTER, FMTl, ITER, FOR  I ,=  (l.B,I,UR)  PO  R(I-1))| 
FOR  J ,a  M STEP  -1  UNTIL  0 DO 
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WRITE  (PRINTER, FMT2,  S(J)  ,FOR  I |»  (LB,1,UB)  DO  U(I»1,J>)| 

LB  I a UP  ♦ 1 | 

END!  •] 

WRITE (PRINTER, F*T3, TEST) f 
END 

•FOR, IBZX  ITERAT 

SUBROUTINE  ITERAT  ‘| 

IMPLICIT  DOUBLE  PRECISION  (C,0,T,U,V) 

COMHON/COEFO/CO(1) 

COMMON/COEF1/C1  (1) 

C0MM0N/C0EF2/C2  ( 1 ) 

C0MMPN/C0EF3/CS(t) 

COHMON/UNK/U ( 1 ) 

COWMON/PAR/ I*AX,N,M, MOOIT, OMEGA, TEST, TEST2 

IF  (HOOIT  ,EQ.  0 ) TEST  a 0,000 

DM  J s 2,M 

DO  1 I a 2,n 

IJ  • I A IMAX  a (J-l) 

IM1J  a IJ  - l 
IP1J  a IJ  ♦ l 
IJM1  s IJ  - IMAX 
UP!  d IJ  ♦ IMAX 
UOLO  a U(IJ) 

UNEW  a -(f.J(IJ)  * Cl  (IJ)*U(IP1J)  * C2(I.I)  *U(IJP1) 

1 ♦ C1(IM1J)*U(IM1J5  + C2(IJM1)*U(IJM1))  / CO(IJ) 

VINT  a (1,000  - OMEGA)  * UOLO  ♦ OMEGA  * UNEW 
U (IJ)  a OMAX1  (VINT, 0,000) 

IF  (HODIT  , ME,  0)  GO  TO  1 
VABS  a DABS  (U(IJ)  - UPLD) 

TEST  a DMAX1  CTEST,VAHS) 

TEST2  a TESTS  ♦ VABS**2 
I CONTINUE 
RETURN 
END 


»XQT 
0 6 
PXOT 
B 12 
PXQT 
16  20 
• XQT 
32 


a.st&o 

72.0UO 

ae.ORRO 

12, OERO 

l,CRR«6 

O.BR&O 

72.0RR0 

4 8 , 0 R R 0 

1 2 , OR  EO 

1 ,0Rtr6 

0.8R&0 

72.0R&0 

48,0RR0 

1 2, ORRP 

1 , 0 E R*6 

R.BftSO 

72.0RR0 

OP, onto 

12.0U0 

1 , 0RR*6 
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Appendix  C: 


Tho  strong  form  of  Green 'r.  theorem 


In  deriving  variational  inequalities  it  is  often  necessary  to  integrate  by  parts 
(i.e.  usq  Green's  theorem).  The  weak  form  of  Green's  theorem  (for  domains  bounded  by 
piecewise  smooth  boundaries)  is  well  known  (Kellogn  [1?53,  p.  84)).  The  strong  form  of 
Green's  theorem  (for  domains  bounded  by  Jordan  curves)  is  less  well  known.  Since  it  is 
desirable  to  make  the  weakest  possible  a-priori  assumptions  about  the  free  boundary,  the 
strong  form  of  Green's  theorem  is  of  value  and  we  therefore  describe  it  here. 

We  recall  that  a closed  Jordan  curve  is  a mapping 

s ■*  z(s)  = (x(s)  , y (s) ) 

of  the  interval  [0,1]  into  the  xy-plane  such  that  z(0)  = z(l)  and  z(s^)  = z(s2>  iff 
either  s^  = s2  or  s = 0 and  s2  = 1 . The  curve  is  rectifi able  if  the  mapping  is  of 
bounded  variation;  that  is,  there  exists  a constant  L such  that  for  all  subdivisions 


0=sn<s1<-'-<E  = 1 > 

0 1 n 


n-1  ..  - ? 

l7-<Ei+1>  - 2<Vl  =.L  l<x(si+i)  - x(s.))  +(y(s.+1)  - y(s.)) 


n-1 

-l 

i=0 
< L . 


2,1/2 


If  J is  a closed  rectifiable  Jordan  curve  then  x(s)  and  y(s)  are  of  bounded  variation 
so  that  the  Riemann-Stielt jes  integrals 

/ f dx  and  / f dy 
J J 

are  defined  for  every  continuous  function  f . 

Now  let  Q be  a bounded  domain  in  the  xy-plane  with  boundary  3f)  which  is  a closed 
rectifiable  Jordan  curve.  Green's  formula  takes  the  form 

//  (|r  “ >d*dy  = / Mdy  + h’dx  , 

n y mi 


where  the  integral  over  3Q  is  taken  in  the  positive  direction  around  3fl  . If  x(s) 
and  y(s)  are  absolutely  continuous  then  the  formula  takes  tho  more  familiar  form 


//  (f“  - ) dxdy  = / (My  + Nx)ds  , 

(t  dX  3y  3fi 

= / (Mil  _ Nn  ) d t. 

Ml  X y 
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where  i denotes  length  and  (n  ,n  ) is  the  unit  outward  normal  which  is  defined  a.c. 

x y 

To  establish  conditions  under  which  Green's  formula  holds,  it  suffices  to  considcr 
the  case  when  one  of  the  functions,  N say,  is  zero. 


Theorem  C.l 


Let  M be  continuous  on  ft  . Assume  that 
(i)  exists  a.e.  and  is  Lebesguo  summable 

Riemann  integrable.) 


3M 

3x 


exists  everywhere  and  is 


(ii)  //  dxdy  = / Mdy 

R SR 

for  every  rectangle  R:  a<x<A,  b<y<B  contained  in  ft,  where  the  integral  over 
R is  taken  in  the  sense  of  Lebesgue  (Riemann) . 

Then  Green's  formula  holds: 


II 

ft 


f)M 

•r-  dxdy  / Mdy 
3ft 


the  integrals  being  in  the  sense  of  Lebesgue  (Riemann) . □ 

Proofs  of  the  above  theorem  are  given  by  Verblunsky  (1949)  for  the  Lebesgue  case 
and  Potts  [1951]  for  the  Riemann  case;  these  authors  give  references  to  earlier  proofs. 

If  condition  (i)  of  the  theorem  is  assumed  then  sufficient  conditions  for-  condition 
(ii)  to  hold  are: 

(a)  M is  absolutely  continuous  in  a ^ x A for  almost  all  y in  (b,B) 

(Lebesgue  case) . 

(3)  The  integral  / Mdy  exists,  (Riemann  case) . 

3r 

In  particular,  if  u is  absolutely  continuous  on  ft  and  continuous  on  ft  , and 
and  if  and  belong  to  L2  (ft)  , then 

II  (!x  " l^)dxdy  = / udy  + udx  . 

ft  dX  dy  3ft 


1 


Bibliography 


Each  item  in  the  bibliography  is  followed  by  a list  of  the  pages  on  which  it  is 
quoted. 


ADAMS,  R.  A.:  Sobolev  Spaces.  Now  York:  Academic  Press,  1975.  [5;17;37;40] 

BABBITT,  H.  E.  and  CALDWELL,  D.  H. : The  free  surface  around,  and  interference  between, 
gravity  wells.  Bulletin  Series  No.  374,  Engineering  Experiment  Station- 
University  of  Illinois,  Urbana,  111.,  1948.  [34] 

BAIOCCHT , C. : Sur  un  probleme  a frontiere  libre  traduisant  le  filtrage  do  liquides  a 

travers  des  milieux  poreux.  Coinptes.  Rendus  Acad.  Sci.  Paris  A273,  1215-1217  (1971).  13] 
BAIOCCHI,  C. : Problemos  a frontiere  libre  et  inequations  variationnelles.  Comptes 
Rendus  Acad.  Sci.  Paris  283A,  29-32  (1976).  110] 

CO  A 

BAIOCCHI,  C. : Estimation  d'erreur  dans  L pour  les  inequations  a obstacle.  Publ  l ic  azioni 
No.  79,  Laboratorio  di  Analisi  Numorica  del  C.N.R.,  Universita  di  Pavia  (1976a). 

[25] 

BAIOCCHI,  C.  , BREZZI,  F.  and  COIIINCIOLI,  V.  : Free  boundary  problems  in  fluid  flow  through 
porous  media.  In  Proceedings  Second  International  Symposium  on  Finite  Element 
Methods  in  Flow  Problems,  Santa  Margherita,  Italy,  1976.  [ 3 J 
BAIOCCHI,  C.,  COMINCIOLI,  V.,  GUERR1 , L. , and  VOLPI,  G. : Free  boundary  problems  in  the 
theory  of  fluid  flow  through  porous  media:  a numerical  approach.  Calcolo  10, 

1-85  (1973).  [4] 

BAIOCCHI,  C.,  COMINCIOLI,  V.,  MACENES,  E. , and  P0Z7.I,  G.  A.:  Free  boundary  problems  in 
the  theory  of  fluid  flow  through  porous  media:  existence  and  uniqueness  theorems. 

Annali  di  Matematica  4 (97) , 1-82  (1973).  ]4;15;16) 

BAIOCCHI,  c.  and  FRIEDMAN,  A.:  A filtration  problem  in  a porous  medium  with  variable 
permeability.  Annali  di  Mat.,  to  appear.  [16:18] 

BEAR,  J. : Dynamics  of  Fluids  in  roious  Media.  New  York:  American  Elsevier,  1972. 

U;4] 


BENCI,  V.:  Su  un  problema  di  filtrazione  in  un  mezzo  poroso  non  omogeneo.  Rend.  Acad. 

Naz.  Lincei  (8)  54,  10-15  0973).  [4;10;15] 

BENCI,  V. : On  a filtration  problem  through  a porous  medium.  Annali  di  Matem.  (4)  100, 
191-209.  (1974).  |4 ; 10; 15; 16; 17 ; 18] 


-49- 


A 


BOULTON,  N.  S.:  The  flow  pattern  near  a gravity  well  in  a uniform  water  bearing  medium. 

J.  Instn.  Civil  Engrs.  3_6,  534-550  (1951)  . (34) 

BREZIS , H.:  Seuil  de  regularite  pour  certains  problemes  unilateraux.  Comptes  Uendus 
Acad.  Sci.  Paris  273A,  35-37  (1971).  [18] 

BREZIS,  H.  and  KINDERLEHRER,  D. : The  smoothness  of  solutions  to  non-linear  variational 
inequalities.  Indiana  Univ.  Math.  J.  23,  831-844  (1973/74).  118) 

* 

BREZIS,  H.  and  STAMPACCHIA,  G.  : Sur  la  regularite  de  la  solution  d' inequations  elliptiques. 

Bull.  Soc.  Math.  France  96,  153-180  (1968).  [18] 

BREZZI,  F. , HAGER,  W.  H. , and  RAVIART,  P.  A.:  Error  estimates  for  the  finite  element 

solution  of  variational  inequalities  Part  I:  Primal  theory.  To  appear.  [18; 21; 25) 
BREZZI,  F.  andSACCHI,  G.  : A finite  element  approximation  of  a variational  inequality 
related  to  hydraulics.  Pubblicazion.i  No.  97,  Laboratorio  di  Analisi  Numerica  del 
C.N.R.,  Universita  di  Pavia,  to  appear.  [21; 25;36; 39] 

CAFFARELLI,  L.  A.:  The  smoothness  of  the  free  surface  on  a filtration  problem.  To 
appear.  [18] 

CIARLET,  P.  G.  and  RAVIART,  P.  A.:  General  Lagrange  and  Hermite  interpolation  in  Rn 

with  applications  to  finite  element  methods.  Arch.  Rational  Mech.  Anal.  46^,  177- 
199  (1972).  [23 ; 40) 

COMINCIOLI,  V.:  Metodi  di  rilassamento  per  la  minimizzazione  in  uno  spazio  prodotto. 
Publicazioni  No.  20,  laboratorio  di  Analisi  Numerica  del  C.N.R.,  Universita  di 
Pavia,  (1971)  . [31] 

COMINCIOLI,  V.:  A theoretical  and  numerical  approach  to  some  free  boundary  problems. 

Annali  di  Matematica  (4)  100,  211-238  (1974).  [34] 

COMINCIOLI,  V.:  On  some  oblique  derivative  problems  arising  in  the  fluid  flow  in  porous 
media.  A theoretical  and  numerical  approach.  Appl.  Math,  and  Optim.  JL,  313-336 
(1974a).  [34] 

COMINCIOLI,  V.:  A comparison  of  algorithms  for  some  free  boundary  problems.  Pubbl icazioni 

No.  79,  Laboratorio  di  Analisi  Numerica  del  C.N.R.,  Universita  di  Pavia,  (1974b). 

[34] 


— 


- 


COMINCIOLI,  V.:  Sur  1 ' approximation  numerique  der,  problcmes  a frontiere  libro  lies  a la 

filtration  dans  les  materiaux  poreux.  Publicazione  No.  105,  Laboratorio  di  Analisi 
Numerica  del  C.N.R.,  Universita  di  Pavia  (1975).  (34) 

COMINCIOLI,  V.,  GUERR1,  L.  , and  VOLPI,  G.  j Analisi  numerica  di  un  problema  di  frontiera 
libera  connosso  col  moto  di  un  fluido  attravorso  un  mezzo  poroso.  Pubbl icazioni 
No.  17,  Laboratorio  di  Analisi  Numerica  del  C.N.R. , Universita  di  Tavia  (1971). 

134) 

COTTLE,  R.  W. : Complementarity  and  variational  problems.  Technical  Report  No.  SOL 
74-6,  Department  of  Operations  Research,  Stanford  University  (1974).  ( 30 ; 34 ] 
COTTLE,  R.  W. : Computational  experience  with  large-scale  linear  complementarity 

problems.  Technical  Report  No.  SOL  74-13,  Department  of  Operations  Research, 
Stanford  University  (1974a) . [30) 

COTTLE,  R.  W.,  GOLUB,  G.  H.,  and  SACHER,  R.  S.:  On  the  solution  of  large,  structured 

linear  complementarity  problems:  III.  Technical  Report  No.  SOL  74-7,  Department  of 
Operations  Research,  Stanford  University  (1974).  [30] 

COTTLE,  R.  W.  and  PANG,  I.-S.:  On  solving  linear  complementarity  problems  as  linear 
programs.  Technical  Report  No.  SOL  76-5,  Department  of  Operations  Research, 
Stanford  University  (1976).  [30] 

COTTLE,  R.  VJ.  and  SACHER,  R.  S.:  On  the  solution  of  large,  structured  linear  comple- 
mentarity problems:  I.  Technical  Report  No.  SOL  73-4,  Department  of  Operations 
Research,  Stanford  University  (1973).  (30] 

COURANT,  R.  and  HILBERT,  D. : Methods  of  Mathematical  Physics,  vol.  2.  New  York: 
Interscience,  1962.  [13] 

CRYER,  C.  W. : The  solution  of  a quadratic  programming  problem  using  systematic  over- 
relaxation. SIAM  J.  Control  9_,  385-392  (1971).  [31] 

CRYER,  C.  W. : A survey  of  steady-state  porous  flow  free  boundary  problems.  Technical 
Summary  Report  No.  1657,  Mathematics  Research  Center,  University  of  Wisconsin, 
Madison,  Wisconsin,  1976.  [1,3,4] 

DURAND,  J-F. : Resolution  numerique  de  problemes  aux  limites  sous-harmoniques . Thesis, 
Univorsite  de  Montpellier,  1968/69.  [30] 


-51- 


DURAND,  J-F. : L' algorithme  de  Ganns-Seidel  applique  a un  probleme  unilateral  non 

symetrique.  Revue  Francaise  d 1 Automat ique , Informatique , ct  Recherche  Operation- 
nelle  6,  No.  R-2,  23-30  (1972).  [30) 

ECKHARDT,  U.  : Quadratic  programming  by  successive  overrclaxation.  Berichte  Nr. 

Jul-1064-MA,  Zentralinstitut  fur  angewandte  Mathematik,  Kernforschungsanlage  Julich, 
Germany,  1974.  [31) 

ELLIOTT,  C.  M. : Some  applications  of  the  finite  element  method  in  numerical  analysis. 

D.  Phil.  Thesis,  Oxford  University,  September  1970.  [4] 

FALK,  R.  S. : Error  estimates  for  the  approximation  of  a class  of  variational  inequalities. 
Math.  Comp.  28,  963-971  (1974).  (25) 

FRANCE,  P.  W. , PAREKH,  C.  J. , PETERS,  J.  C. , and  TAYLOR,  C. : Numerical  analysis  of  free 
surface  seepage  problems.  J.  Irrigation  Drainage  Division,  Proc.  Amer.  Soc.  Civil 
Engrs.  97  No.  IR1,  165-179  (1971).  [35] 

FRIDMAN,  V.  M.  and  CHERNINA,  S.  V. : An  iteration  process  for  the  solution  of  the  finite 
dimensional  contact  problem.  U.S.S.R.  Comput.  Math,  and  Math.  Physics  1_,  210-214 
(1967;.  =Zh.  Vych.  Mat.  mat.  Fiz  7,  160-163  (1967).  [30] 

FRIEDMAN,  A.:  Personal  communication,  May  1977.  [16] 

GLOWINSKI,  R.  : La  methode  de  relaxation.  Rendiconti  di  Matematica  1^.<  Universita  Roma, 
(1971).  [30] 

GLOWINSKI,  G. : Sur  la  minimization,  par  surrolaxation  avec  projection  de  fonctionnelles 
quadratiquos  dans  la  cspaces  de  Hilbert.  Comptes  Rendus  Acad.  Sci.  Paris  276 , 
1421-1423  (1973).  [30] 

GLOWINSKI,  R. : Introduction  to  the  approximation  of  elliptic  variational  inequalities. 
Report  of  the  Laboratoirc  Analyse  Numerique,  Univcrnite  Paris  VI,  1976.  [25] 

GLOWINSKI,  R. , LIONS,  J.  L.  and  TREMOLIERES,  R. : Analyse  Numerique  des  Inequations 
Var iationnel les ■ Paris:  Dunod,  1976.  [25] 

HAGER,  W.  W. : Rates  of  convergence  for  discrete  approximations  to  unconstrained  control 
problems.  SIAM  J.  Numer.  Anal.  1_3,  449-472  [1976].  [25] 

HAGER,  W.  W.  and  STRANG,  G. : Free  boundaries  and  finite  elements  in  one  dimension. 

Math.  Comp.  29,  1020-1031  (1975).  (25) 

HALL,  H.  P. : An  investigation  of  steady  flow  towards  a gravity  well.  La  Houille  Blanche 
10,  8-35  (1955) . [33; 34) 


52- 


HANTUSH,  M.  S.:  Hydraulics  of  walls.  Advances  of  Hydroscience  1,  281-432.  (1964). 

Ui4] 

HERBERT, ' R. : Time  variant  groundwater  flow  by  resistance  network  analogues.  J.  Hydrology 
6,  237-264  (1968).  (39] 

HILDRETH,  C.  G.  : A quadratic  programming  procedure.  Naval  Res.  Legist.  Quart.  4_,  79-85 
[1957).  [30] 

JENSEN,  R.  : Structure  of  the  non-monotone  free  boundaries  in  a filtration  problem.  To 
appear.  [18] 

KASHEF,  A.:  Numerical  solutions  of  steady-state  and  transient  flow  problems,  artesian 
and  water  well  problems.  Fh.D.  thesis,  Purdue  University  (1953).  [34] 

KASHEF,  A.  I,,  T0U1.0UKIAN,  Y.  S.  , and  FADUM,  R.  E.  : Numerical  solution  of  steady-state 
and  transient  flow  problems  - artesian  and  water-table  wells.  Purdue  University 
Engineering  Experiment  Station  Bull.  117,  Lafayette,  Indiana  (1952).  [34] 

KEALY , c.  D.  and  BUSCH,  R.  A.:  Determining  seepage  characteristics  of  mill-tailings 

dams  by  the  finite-element  method.  Rep.  Invest.  RI  7477,  U.  S.  Bureau  of  Mines, 
Washington,  D.  C.,  January  1971.  [35] 

KELLOGG,  O.  D.  : Found. -1: ions  of  Potential  Theory.  New  York:  Dover,  1953.  [47] 

KIRKHAM,  D. : Exact  theory  for  the  shape  of  the  free  water  surface  about  a well  in  a 
semiconfincd  aquifer.  J.  Geophysical  Res.  69,  2537-2549  (1964).  [34] 

LEWY,  H.  and  STAMPAGCH1A,  G. : On  the  regularity  of  the  solution  of  a variational  in- 
equality. Communications  Pure.  Appl . Math.  72-'  153-188  (1969).  [18] 

LIONS,  J.  L. : Optimal  Contiol  of  Systems  Governed  by  Partial  Differential  Equations. 
Berlin:  Springer,  1971.  [20] 

LUONG,  N. : Sur  la  metbode  de  sur-relaxation,  dans  le  cas  dcs  problemes  avec  contraintc: 
un  resultat  de  convergence  asymptotiquo.  Rev.  Francaise  Automat.  Informat. 

Recherche  Operationnele  2.  N°*  R-2,  107-113  (1973).  [31] 

MANGASAKI AN , 0.  L. : Solution  of  symmetric  linear  complementarity  problems  by  iterative 
methods.  Technical  Report  No.  275,  Computer  Sciences  Dept.,  University  of 
Wisconsin,  Madison,  Wisconsin,  1976.  [31] 

MARTINET,  B. : Convergence  do  certaines  methodes  de  relaxation  on  progr ammati on  convexe. 
Comptes  Rendus  Acad.  Sci.  Paris  265A,  210-212  (1967).  [30] 


-53- 


MARTINET,  B.  and  AUSLENDER,  A.:  Methodes  de  decomposition  pour  la  minimisation  d'une 
fonction  sur  une  espacc  produit.  SIAM  J.  Control,  12_,  635-642  (1974).  (311 

MAUERSBERGER,  P. : Ein  Variationspri nzip  fur  die  stationare  Grundwasserstromung  in  der 

Umgebung  eines  vollkommenen  Brunnens.  Gerlands  Beitrage  zur  Geopliysik  74,  343-350 
(1965)  . [4] 

MAUERSBERGER,  P.:  Bemerkungen  zur  Theorie  des  stationaren  Brunnens.  Gerlands  Beitrage 
zur  Geophysik  74,  516-526  (1965a) . [4] 

MAUERSBERGER,  P.:  Herleitung  einer  ersten  Naherung  fur  Grundwasserspiegel  und  Stromungsfeld 
eines  stationaren  Brunnens  mit  Hilfe  der  Trefftzschen  Variationsmethode.  Acta 
Hidrophysica  11_,  171-179  (1967).  [35] 

MAUERSBERGER,  P. : The  use  of  variational  methods  and  of  error  distribution  principles  in 
groundwater  hydraulics.  Bull.  Internat.  Assoc.  Scientific  Hydrology  ]_3,  169-178 
(1968).  [35] 

MAUERSBERGER,  P.:  Eine  hydrologische  Anwendung  der  Trefftzschen  Variationsmethode. 

Gerlands  Beitrage  zur  Geophysik  22'  235-250  (1968a).  [35] 

MAUERSBERGER,  P.  : Die  Methode  der  Randkollokation  in  der  Theorie  stationarer  Grundwasser-- 
bewegungen.  Gerlands  Beitrage  zur  Geophysik  22'  331-336  (1968b).  135] 

MAUERSBERGER,  P.:  Fehlerabgleischsprinzipien  in  der  Theorie  stationarer  dreidimensionaler 
Grundwasserbewegungen  mit  freier  Oberflache.  Gerlands  Beitrage  zur  Geophysik  77 , 
363-374  (1968c).  [35] 

MAUERSBERGER,  P. : Die  Ergiebigkeit  eines  stationaren  Brunnens  im  inhomogenen,  orthotropen 

Grundwassertrager . Acta  Hydrophysics-  14,  143-163  (1969).  [4, -16] 

MERZLJAKOV,  Ju.  I.:  On  a relaxation  method  for  solving  systems  of  linear  inequalities. 

USSR  Comput.  Math,  and  Math.  Phys.  482-487  (1962).  [30] 

MIELLOU,  J.-C.  : Methodes  de  Jacobi,  Gauss-Scidel,  sur-(sours)  relaxation  par  blocs, 
appliquees  a une  classc  de  problemes  non  lineaircs.  Comptes  Rendus  Acad.  Sci. 

Paris  A273 , 1257-1260  (1971).  [31] 

MIELLOU,  J.-C.:  Methodes  de  Gauss-Scidel , et  sur  relaxation  par  blocs,  dans  le  cas  de 

problemes  non  lineaircs,  application  au  controle  optimal  de  systemes  gouvernes  par 
des  equations  aux  derivees  partielles:  methode  de  l'etat  adjoint  par  relaxation. 
Expose  au  Colloque  d'Analyse  Numerique  d'Anglct,  1971a.  [31] 


-54- 


MIELLOU,  J.-C.:  Sur  une  variante  de  la  methode  de  relaxation  appliquee  a des  problomcs 
comportant  un  operateur  somme  d'un  operateur  differentiable,  et  d'un  operateur 
maximal  monotone  diagonal.  Comptes  Rendus  Acad.  Sci . Paris  A275,  1107-1110  (1972).  (31) 

MOSCO,  U.  and  STRANG,  G. : One-sided  approximation  and  variational  inequalities.  Bull. 

Amer.  Math.  Soc.  80,  303-312  (1974).  [25] 

MURRAY,  J.  A.:  Relaxation  methods  applied  to  seepage  flow  problems  in  earth  dams  and 
drainage  wells.  J.  Instn.  Engrs.  (India)  £1,  No.  4,  149-161  (1960).  l34l 
NEUMAN,  S.  P.  and  WITHERSPOON,  P.  A. : Finite  element  method  of  analyzing  steady  seepage 
with  a free  surface.  Water  Resources  Res.  6,  889-897  (1970).  133; 34 1 
NEUMAN,  S.  P.  and  WITHERSPOON,  P.  A.:  Analysis  of  nonsteady  flow  with  a free  surface 
using  the  finite  element  method.  Water  Resources  Res.  T_,  611-623  (1971).  [33;  34 ) 
POLUBARINOVA-KOCHINA,  P.  Ya.:  Theory  of  Ground  Water  Movement.  Princeton:  Princeton 
University  Press,  1962.  (4) 

POTTS,  D.  H. : A note  on  Green's  theorem.  J.  London  Math.  Soc.  26^  302-304  (1951).  (48) 

SCHMIDT,  II.:  Uber  eine  Anwendung  der  Relaxiationsmethode  zur  Behand)ung  von  Grund- 
wasserstromungon.  Diss.  T.  H.  Wien,  1956.  (34) 

STAMPACCHIA,  G. : Formes  bilinaaires  coercitivcs  sur  les  ensembles  convexes.  Comptes 
Rendus  Acad.  Sci.  Paris  258,  4413-4416  (1964).  (17] 

TAYLOR,  G.  S.  and  LUTHIN,  J.  N.:  Computer  methods  for  transient  analysis  of  water-table 
aquifers.  Water  Resources  Res.  !>,  144-152  (1969).  (33:34] 

TAYLOR,  R.  L. : Axisymmetric  and  plane  flow  in  porous  media.  Technical  Report,  University 
of  California  at  Berkeley,  1966.  (34J 

TAYLOR,  R.  L.  and  BROWN,  C.  B. : Darcy  flow  solutions  with  a free  surface.  J.  Hydraulics 
Division,  Proc.  Amer.  Soc.  Civil  Engrs.  93,  No.  1IY2,  25-33  (1967).  (351 

VARGA,  R.  S. : Matrix  Iterative  Analysis.  Englewood  Cliffs:  Prentice-Hall  (1962).  I3H 
VERBLUNSKY,  S. : On  Green’s  formula.  J.  London  Math.  Soc.  2 £,  146-148  (1949).  (481 

YANG,  S.  T. ; Seepage  towards  a well  analyzed  by  the  relaxation  method.  Ph.D.  thesis, 

Harvard  University,  (1949).  (34) 

YOUNGS,  E.  G. : Steady  state  flow  around  wells  in  aquifers  with  hydraulic  conductivity 

varying  with  depth.  Water  Resources  Res.  1366-1368  (1971)  . [Errata:  According 
to  the  author  the  figures  in  Table  1 are  slightly  in  error].  [4;16] 


-55- 


SECURITY  CLASSIFICATION  of  This  pace  (Whan  Data  Fntarad) 


REPORT  DOCUMENTATION  PAGE 


report  NUMBER 


2.  GOVT  ACCESSION  NO. 


READ  INSTRUCTIONS 
BEFORE  COMPLETING  FORM 

3.  RECIPIENT'S  CAT  ALOG  NUMBER 


4.  TITLE  (and  Subtitle) 

THE  NUMERICAL  SOLUTION  OF  AXISYMMETRIC  FREE 
BOUNDARY  POROUS  FLOW  WELL  PROBLEMS  USING 
VARIATIONAL  INEQUALITIES 

7.  author^; 

Colin  W.  Cryer  and  Hans  Fetter 

• PERFORMING  ORGANIZATION  NAME  AND  AOClRCSS 

Mathematics  Research  Center^  University  of 

610  Walnut  Street  . Wisconsin 

Madison.  Wisconsin  53706 

II.  CONTROLLING  OFFICE  NAME  AND  ADDRESS 

See  Item  18  below. 


5.  type  of  report  a perioo  covered 

Summary  Report  - no  specific 
reporting  period 

6.  PERFORMING  ORG.  RCPORT  NUMBER 
V CONTRACT  OR  GRANT  NUMBLRf a) 

/ 

DAAG29-7  5-C-0024 

DCR75-03838 

io.  phograh“ellmi:nt  projlct,  t ask 

AKEA  A WORK  UNIT  NUMBERS 

Work  Unit  Number  7 - 
Numerical  Analysis 

12.  REPORT  DATE 

June  1977 

13.  NUMBLR  OF  PAGES  , 


MONITORING  AGLNCY  KAMI:  A AC'DKLSSfU  dilletatii  from  Controlling  Ottico)  IS.  SECURITY  CLASS,  (e/  thi a report) 

UNCLASSIFIED 

~Wm~  t'ifc  CL  ASSl  Fj  C ATIOli/ DOWNGRAOIN  G*“ 
SCHEDULE 

16.  DISTRIBUTION  STA»  EMEIlT  (ol  ihla~ oil) 

I Approved  for  public  release:  distribution  unlimited. 


I 17.  DISTRIBUTION  ST  AT  r MEN T (ol  lha  ebatcbcl  antarad  fn  Bloc*  20,  It  d.ttarant  l:o w Import) 


l»  SUPPLtMENT ARY  NOTES 

U.SJ  Army  Research  Office 

r.O.  Box  122X1  and  Nations 

Research  Triangle  Park  Washing' 

North  Carolina  27709 

19.  KEY  WORDS  (Continue  on  ravotaa  aioo  II  nae eaaa:r  wnd  Idantlty  Ly  block  ntrmber) 


National  Science  Foundation 
Washington,  D.  C.  20550 


Porous  flow.  Free  boundary  problems.  Unconfined  flow,  Numerical  methods, 
Variational  inequalities.  Existence,  Uniqueness,  Error  estimates.  Finite 
, elements  , a " f 

^ a c'~<e 

' 2 

m/aBSTRACT  (Contlnuo  on  '.•*,*■  alda  II  nmcoobmay  mild  Idantlty  by  block  nuaib.r) 

The  free  boundary  problem  for  a fully  penetrating  well  of  radius  r and  filled 
with  water  to  a depth  h,  in  a layer  of  soil  of  depth  U,  ^radius  R and  per- 
meability k(x,y)  can  be  formulated  as  follows:  Find  <f c C |r,R]  and  u e C (fi) 

n C(H)  such  that  (*kux)x  + (xkuy)y  = 0 in  0,  u(R,y)  - H for  0 < y < H , 

u (x,0)  = 0 for  r < x < R,  u(r,y)  = h for  0 <y  < h,  u(r,y)  = y for  h<  y <<?(r)  , 

u (x,i P (x) ) = 0 and  u(x,^  (x) ) = V (x)  for  r < x < R,  where  £1  = 0x«y)  : t < x < R, 
0*<y<^(x)J  . The  results  of  Benci  {Annali  di  Mat.  100(1974),  191-209)  are  use 


0 , J*H  7j  1473  ED'T,C,N  or  ’ MOV  “ *s, obsolete  UNCLASSIFIED 

SECURITY  CLASSIFICATION  OF’thIS  PALE  flHi.IT  D.l.  Knl.r.J) 

problem  is  approximated  using  piecewise  linear  finite  elements  and  0(h)  convergence 
of  the  approximate  solutions  is  proved  using  recent  results  due  to  Brezzi,  Hager, 
and  Kaviart^.  ) 


